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We consider the two dimensional (2D) classical lattice Coulomb gas as a model for magnetic field 
induced vortices in 2D superconducting networks. Two different dynamical rules are introduced 
to investigate driven diffusive steady states far from equilibrium as a function of temperature and 
driving force. The resulting steady states differ dramatically depending on which dynamical rule is 
used. We show that the commonly used driven diffusive Metropolis Monte Carlo dynamics contains 
unphysical intrinsic randomness that destroys the spatial ordering present in equilibrium (the vortex 
lattice) over most of the driven phase diagram. A continuous time Monte Carlo (CTMC) is then 
developed, which results in spatially ordered driven states at low temperature in finite sized systems. 
We show that CTMC is the natural discretization of continuum Langevin dynamics, and argue that 
it gives the correct physical behavior when the discrete grid represents the minima of a periodic 
potential. We use detailed finite size scaling methods to analyze the spatial structure of the steady 
states. We find that finite size effects can be subtle and that very long simulation times can be 
needed to arrive at the correct steady state. For particles moving on a triangular grid, we find that 
the ordered moving state is a transversely pinned smectic that becomes unstable to an anisotropic 
liquid on sufficiently large length scales. For particles moving on a square grid, the moving state is 
a similar smectic at large drives, but we find evidence for a possible moving solid at lower drives. 
We find that the driven liquid on the square grid has long range hexatic order, and we explain this 
as a specifically non-equilibrium effect. We show that, in the liquid, fluctuations about the average 
center of mass motion are diffusive in both the transverse and longitudinal directions. 
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I. INTRODUCTION 

While the theory of phase transitions of systems in 
thermodynamic equilibrium is a well established and ma- 
ture area of statistical physics, much less is established 
about analogous critical behavior in driven steady states 
far from equilibrium. As has been the case in the study of 
equilibrium phase transitions, the use of lattice models, 
in which the degrees of freedom are constrained to sit on 
the sites of a discrete periodic grid, has led to analytical 
simplifications and greater accuracy in numerical simu- 
lations for investigating such steady states, as compared 
to corresponding continuum modelsi. One advantage of 
a lattice gas model for numerical simulations of driven 
interacting many-body systems is that particles hop in 
discrete jumps. If a particle sits in a local potential min- 
imum, the lattice gas dynamics can allow the particle 
to hop over the energy barrier out of the minimum in a 
single move. In contrast, in continuum simulations like 
molecular dynamics, considerable simulation time can be 
wasted at low temperatures waiting for a thermal excita- 
tion that will excite the particle over the energy barrier. 
The lattice gas method can therefore hope to simulate 
out to much longer effective times, and focus on effects 
due to many-body interactions rather than single body 
potentials. 

One of the first, and still one of the most commonly 
used, numerical methods to simulate driven steady states 
of a lattice gas is the driven diffusive Monte Carlo 
method. This method, introduced by Katz, Leibowitz 



and Spohn 2 , extends familiar equilibrium Monte Carlo 
methods to the case of driven non-equilibrium states. 
The key idea of this method is to include the work done 
by the driving force on a moved particle, in addition to 
the change in interaction energy, when computing the 
energy difference to use in the Monte Carlo test for mak- 
ing moves. Such a term biases motion in the direction of 
the driving force, and, with the use of periodic boundary 
conditions, results in a steady state with a finite parti- 
cle current. This algorithm, which satisfies local detailed 
balance for individual particle moves, seeks to model dif- 
fusively moving particles in the limit where motion is 
dominated by thermal activation over energy barriers, 
rather than microscopic dynamics. The hope is that the 
main qualitative features of the driven steady states, and 
possible phase transitions between them, will be indepen- 
dent of the details of the microscopic dynamics, and so 
will be captured by this algorithm. 

However, unlike equilibrium simulations, where any 
dynamics that satisfies detailed balance is sufficient to 
generate the correct equilibrium Gibbs ensemble and so 
equilibrium averages are in principal independent of the 
microscopic dynamics, there is no such guarantee for 
non-equilibrium states. Even for sets of dynamics that 
would appear to lie within the same dynamic "universal- 
ity class"— from the viewpoint of symmetry and conserved 
quantities, averages in driven steady states far from equi- 
librium may conceivably be qualitatively, not just quan- 
titatively, different. 

In this work we test this notion explicitly by consider- 
ing two different versions of driven diffusive Monte Carlo 
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dynamics, both intended to model the overdamped diffu- 
sive limit. We consider first (i) driven diffusive Metropolis 
Monte Carlo dynamics-*- (DDMMC), where the standard 
Metropolis method is used to select attempted excita- 
tions and decide whether or not to accept them. We then 
consider (ii) driven diffusive continuous time Monte Carlo 
dynamics (CTMC), where the continuous time Monte 
Carlo method^ is used to make a rejectionless dynam- 
ics. We believe that our work is the first application of 
continuous time Monte Carlo in the context of driven dif- 
fusive problems. We apply these methods to the problem 
of the driven two dimensional (2D) classical one compo- 
nent lattice Coulomb gas, which serves as a model for log- 
arithmically interacting, magnetic field induced, vortices 
in periodic 2D superconducting networks 6 . We consider 
both the cases of a triangular and a square grid of sites. 
This model is of interest because it allows one to consider 
the effect of a uniform driving force on a system which 
has spatially ordered states in equilibrium (the vortex 
lattice), in contrast to simpler nearest neighboring Ising- 
like lattice gas models^, which in general have no such 
periodic spatial order. 

We find that our two dynamics result in dramatically 
different driven steady states, when the system is acted 
upon by a uniform applied force F. We find that over 
most of the T — F phase diagram, the DDMMC method 
results in a spatially disordered moving steady state with 
a very short translational correlation length. We ar- 
gue that this behavior is due to intrinsic randomness in 
the DDMMC algorithm, that is sufficient to disorder the 
moving system even at T = 0. In contrast, we find that 
the CTMC method, at least for finite size systems, can 
result in spatially ordered moving steady states, as well as 
orientationally ordered moving liquids. We demonstrate 
that the CTMC method is the correct discretization of 
diffusive Langevin dynamics in a certain limit, and argue 
that it more generally describes motion when the dis- 
crete grid is thought of as representing the minima of a 
one body periodic potential, and the energy barriers of 
this potential are large compared to the energy change of 
hopping between minima. Thus we believe that CTMC 
is not only a more interesting dynamics, but also a more 
physically correct one. For CTMC dynamics, we carry 
out detailed finite size scaling analyses of our ordered 
steady states, and show that there can be subtle finite 
size effects due to diverging correlation lengths at low 
temperatures. We also show that exceedingly long sim- 
ulations are needed, in some cases, in order to arrive at 
the correct steady state distribution. 

The remainder of this paper is organized as follows. 
In section ITT1 we define in detail our Coloumb gas model 
and our two lattice gas dynamics. We discuss qualita- 
tive behaviors to be expected at low temperatures, and 
define the observables we will measure. In section IIIII 
we present the results of our simulations on a triangular 
grid of sites. We show the phase diagrams of both the 
DDMMC and CTMC for a system of a given finite size, 
and demonstrate the dramatic difference between them. 



We then focus the remainder of our work on CTMC. We 
carry out detailed finite size scaling analyses to study the 
structural order of the moving steady states in both the 
high drive and low drive limits. At low temperature we 
find an ordered moving smectic state, however we argue 
that this state is ultimately unstable to a liquid on suf- 
ficiently large length scales. We also present results for 
dynamic behavior, studying the average velocity of the 
system and the diffusion of the system about its center 
of mass motion. In section IW1 we present our results for 
simulations on a square grid of sites. We study several 
specific points in the phase diagram representative of the 
high drive and low drive limits. Unlike for the triangular 
grid, we find that the structure of the ordered moving 
state appears to have different periodicities at different 
driving forces. We carry out finite size scaling to inves- 
tigate the stability of the ordered states in the large size 
limit. We show that, unlike the liquid state in equilib- 
rium, the liquid driven steady state possesses long range 
hexatic orientational order. In section we discuss our 
results and present our conclusions. 

Some aspects of this work, focusing on the differences 
between DDMMC and CTMC and the structural order 
of driven states on the triangular grid, have previously 
appeared as a letter—. The detailed discussion of the 
phase diagram on the triangular grid, the finite size scal- 
ing analyses, the discussion of dynamical behavior, and 
all results for the square grid, are presented for the first 
time in the current work. 



II. MODEL AND METHODS 

A. Coulomb Gas Model 

Our model is a classical one component lattice 
Coulomb gas of 2D interacting charges, which may be 
taken as a model for interacting vortices in a 2D super- 
conducting network 6 .. The charges are constrained to sit 
on the discrete sites i of a periodic 2D grid. If the ba- 
sis vectors of the grid are {01,62}, we take the grid to 
have finite length L M in direction a M and we take periodic 
boundary conditions in both directions. The Hamilto- 
nian is given by&, 

where the sum is over all pairs of sites i, j of the grid, 
ni G {0,1} is the integer charge on site i, —f is a uni- 
form neutralizing background charge, and G(r) is the 2D 
lattice Coulomb potential which solves, 

A 2 G(r) = -2nS rfi . (2) 

where A 2 is the discrete Laplacian for the grid. Defining 
A 2 appropriate to periodic boundary conditions results 
in a G(r) that satisfies periodic boundary conditions. For 
separations large compared to the grid spacing, but small 
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FIG. 1: (a) Triangular grid of size L\ x L2 with basis vectors 
ai and &2- (b) Reciprocal space to the triangular grid, with 
basis vectors bi and b2. Allowed wavevectors for Fourier 
transforms of real space quantities can be restricted to the 
hexagonal first Brillouin zone shown in (b). 



compared to the grid length (1 <C |r| <C L), one has 
G(r) ~ — In I r I . The condition that the total energy re- 
main finite imposes the charge neutrality condition, 



Nr. 



(3) 



We will consider first the case of a triangular grid of 
sites. Here the basis vectors are, 



a 1 

a 2 



1- \/3\ 

—x H y 

2 2 y 



(4) 



and the sites i of the grid are specified by the position 
vectors, 



{bi, so that they all lie in the hexagonal shaped first 
Brillouin zone of the reciprocal grid. The geometry of 
these allowed wavevectors is illustrated in Fig.^j. 

Defining 03 = a\ — 0,2, the discrete Laplacian for the 
triangular grid is given by, 



A 2 G(r) = c J2 [G(r + a M ) - 2G(r) + G(r - a M )] (8) 



with c an appropriate geometrical constant to give the 
correct continuum limit. Taking the Fourier transform 
of the above, we find that the solution to Eq. (J2J) is given 

by&a, 



G(r) 



7T 

cN 



E 



cos(k ■ di) — cos(k • a 2 ) ~ cos(k • 0,3) 

(9) 

where N = L\Li is the number of sites in the grid, and 
the sum is over all the allowed wavevectors of Eq. J7|). 
The correct value of the geometric constant is c = l/v3. 
However, in order to compare with previous work done 
on this model 6 , we will make the choice c = 2/3. This 
difference amounts only to a rescaling of the temperature. 

We will also discuss the case of a square grid of 
sites. Here the real space basis vectors of the grid are 
{ai, 0,2} = {x, ?/}, the basis vectors of the reciprocal space 
are {bi,b 2 } = {2-Kx,2iry}, and the lattice Coulomb po- 
tential is given by 6 ^, 



ikr 



N ^ 2 



cos(k • Si) — cos(k ■ a 2 ) 



(10) 



m M G {0,1,..., i p - 1} 



(5) 



The geometry of this real space grid is illustrated in 
Fig.DJi. 

The solution to Eq. @ will be given in terms of its 
Fourier transform. The basis vectors of the reciprocal to 
the grid in Fourier transform space are then, 



bi = 2nx - -j=y 

47T A 

b2 = T3 y ' 



(6) 



and the allowed wavevectors consistent with periodic 
boundary conditions are given by, 



fcibi + k 2 b 2 , 



kfj, G {0, 



1 



(T) 



Equivalently, one could translate the above wavevectors 
by an appropriate linear combination of the basis vectors 



B. Lattice Gas Dynamics 

Our goal is to simulate the non-equilibrium steady 
states of the lattice Coulomb gas when driven by a uni- 
form applied force F. For equilibrium simulations, any 
dynamical rule that satisfies detailed balance will con- 
verge to the correct equilibrium ensemble; the details of 
the dynamics may effect the speed of convergence, but 
they are otherwise irrelevant for computing time inde- 
pendent thermodynamic averages. For non- equilibrium 
steady states, however, even time independent averages 
may depend on the details of the microscopic dynamics. 
Here we consider two different microscopic dynamics for 
the simple case of over damped diffusively moving parti- 
cles (the simplest case in the classification scheme of dy- 
namic critical phenomena by Halperin and Hohenbcrg 3 ). 
Both dynamics involve single particle moves only. We 
find that, for finite size systems, the resulting steady 
states for these two dynamics can be qualitatively dif- 
ferent. 
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1. Driven Diffusive Metropolis Monte Carlo (DDMMC) 

The first lattice gas dynamics we consider is the com- 
monly used driven diffusive Metropolis Monte Carlo 
dynamics^ (DDMMC). This algorithm was intro- 
duced as a simple modification of ordinary equilibrium 
Metropolis Monte Carlo. It was intended to model the 
steady states of a driven system in the limit where motion 
is dominated by thermal activation over energy barriers, 
and so presumably is not very sensitive to microscopic 
details. The DDMMC algorithm is defined as follows. 
At each step of the simulation a charge rii = 1 is selected 
at random, and the charge is then moved a test displace- 
ment Ar to a randomly chosen nearest neighbor site. For 
the triangular grid, Ar is chosen with equal probabil- 
ity from the six possible directions ±a M , ^ = 1,2,3. If 
Hold and H. ncw are the interaction energies, Eq. QJ, of 
the system before and after this test move is made, one 
computes the energy difference, 



AU = H„ 



n. 



old 



F • Ar 



(11) 



where the last term is the work done by the applied force 
on the moved charge. This test move is then accepted or 
rejected according to the usual Metropolis criterion, 



accept if 
reject if 



e -AU/T > r 

otherwise , 



(12) 



where r is a random variable uniformly distributed on 
[0, 1). One pass of N c such steps equals one unit of sim- 
ulation time. The term in Eq. I|ll|) proportional to the 
force F biases moves parallel to F and, in conjunction 
with the periodic boundary conditions, will in general set 
up a steady state with a finite current of particles moving 
parallel to F. Time independent averages are computed 
in the usual Monte Carlo way, as a direct average over 
configurations sampled every N pass passes. 



2. Driven Diffusive Continuous Time Monte Carlo 
(CTMC) 

The second dynamics we consider, and the one which is 
used for the main part of our work presented here, we call 
driven diffusive continuous time Monte Carlo (CTMC). 
The algorithm is defined as follows. Starting from a par- 
ticular initial configuration, we denote by {ia) the single 
particle move of a charge n, — 1 on site rj, to its near- 
est neighbor site in direction a. For the triangular grid, 
a G {±ai, ±02, ±03}. For a grid in which each site has 
z nearest neighbors, the total number of such possible 
single particle moves is zN c . For each such move, we 
compute the energy change AUi a according to Eq. (JTTJ, 
and define a probability rate for making this move, 



W e 



-AU ia /2T 



(13) 



where 1/Wo sets the unit of time. Note that the above 
rates, as well as the Metropolis rates of DDMMC set by 



Eq. 112(1. obey a local detailed balance. If s is an initial 
state, and s' is the state reached from s by making the 
single particle move {ia), then, 



W(s -> s') 
W{s> -> s) 



= e -AU ta /T 



(14) 



Although systems out of equilibrium do not in general 
need to satisfy detailed balance, local detailed balance is 
physically reasonable if we wish to regard each charge as 
moving in a local potential determined by its interactions 
with the other charges and with the applied force. 

Having specified the rates of Eq. ((13J1 , we determine 
which move to make by regarding all zN c of the possible 
single particle moves as independent Poisson processes. 
The probability that the next move will be {ia) is then, 



p. = 



and the average time it takes to make this move is. 



At 



1 



1 



(15) 



(16) 



We thus make a move by sampling the probability distri- 
bution Pi a of Eq. (|15fl , and then update the simulation 
clock by the amount At of Eq. (|16fl . Averages of observ- 
ables O are computed as, 

(O) = i J o{t)dt = \Y. ° sAts ' 

where s labels the steps of the simulation, O s is the value 
of O in the configuration at step s, At s = t s+ i — t s is the 
time spent in the configuration at step s according to 
the simulation clock, and r = ^ s At s is the total time 
of the simulation. As in DDMMC, we will refer to N c 
simulation steps as one pass through the system. 

The above algorithm is a straightforward extension of 
equilibrium continuous time Monte Carlo^, but we be- 
lieve that this is its first application in the context of 
driven non-equilibrium steady states. The method was 
first introduced as the "n-fold way" for spin models^. It 
owes its name to the continuous variations in the time 
steps At s , which vary from configuration to configura- 
tion, according to the energy barriers present in each 
configuration. It is a rejectionless algorithm designed 
to speed up equilibration at low temperatures. Rather 
than waste many rejected moves until a rare acceptance 
takes one up and over an energy barrier, the energy bar- 
riers AUia themselves set the time scale for each move, 
which then happens in a single simulation step. Simu- 
lation clock times can vary over orders of magnitude as 
either T or the height of the energy barriers vary. 

In CTMC there are many possible choices for the rates 
Wi a that would satisfy local detailed balance. In Ap- 
pendix A we show that the particular choice of Eq. H13|) , 
with Wq = cDT {D is the diffusion constant), is the 
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natural discretization to a periodic grid of sites of over 
damped continuum Langevin dynamics, and that the 
continuum limit is reached when AU{ a <C T . Our simu- 
lations, however, being generally at low T or large F, are 
mostly in the opposite limit of AUi a > T. To see what 
physical situation this limit corresponds to, consider a 
single particle moving on a one dimensional grid of sites, 
in a driving force F, According to the CTMC algorithm, 
the average distance traveled in one step is, 



(Ax) 



3 F/2T 



-F/2T 



*F/2T 



-F/2T 



while the average time for this step is given by, 



± = W tot = W 



e F/2T + e -F/2T 



leading to an average velocity, 



<Az) 
At 



2W sinh(F/2T) 



(18) 



(19) 



(20) 



At low ratios of F/T the velocity is linear in the ap- 
plied force, (v) ~ WqF/T; at large F/T, the veloc- 
ity grows exponentially, (v) ~ Woe F ' 2T . We can com- 
pare the above result to that of an over damped particle 
moving in a continuum "washboard potential", U(x) = 
— Uo cos(27ra) — Fx, which has been studied by Ambc- 
gaokar and Halperin in the context of a driven Josephson 
junction^. The average velocity that they find agrees ex- 
actly with Eq. H20fl above, if one is in the limit T -C 2Uo 
and F < 2ttUq, and one identifies 

W = 2irU D^T^e~ 2U °{^ T ^ + ~ ts ' m ~ 1 ~ l \ /T , (21) 

where D is the diffusion constant and 7 = F/2ttUq. For 
small 7 the above Wq reduces to a form independent of 
the drive F, 

W ~ 2i:U Q De- 2Uo/T , whenF < 2nU , (22) 

which is the rate for activation over an energy barrier 
of 2Z7o- CTMC thus describes the limit where the grid 
sites represent the minima of a periodic pinning potential, 
and the applied force is weak enough that motion is due 
to thermal activation of particles, one at a time, over 
the barriers of this periodic potential. It is unclear 11 
if CTMC can qualitatively describe the very large drive 
limit, F 3> Uq, where the washboard potential loses its 
local minima parallel to F, and the average velocity again 
becomes proportional to F. For the results reported in 
the following sections we will measure time in units where 
Wq = 1, independent of the temperature T or driving 
force F. 



3. Behavior at Low Temperature 

To get a better feel for the behavior of the above two 
lattice gas dynamic algorithms, we can consider their be- 
havior at low temperature. In the limit T — > 0, the 



DDMMC will reject all moves except those that lower the 
energy, i.e. AUi a < 0. If one starts initially in the F = 
ground state and increases F, all moves will be rejected 
until F reaches a critical value F c equal to the interac- 
tion energy associated with moving one charge forward 
parallel to F. The ordered ground state charge lattice 
will therefore be pinned with (v) = for F < F c , and 
moving with (v) finite for F > F c . 

Next we consider DDMMC at T = with F > F c , 
so that the work done by the force in Eq. (|llfl dominates 
the interaction energy AH. In this case, the DDMMC 
algorithm randomly picks a charge, and then randomly 
picks a direction in which to move it. The move will be 
accepted only if it lowers the energy, i.e. if the charge 
advances in the direction of F. This will happen only 
for a certain fraction p of the possible directions. For a 
triangular grid, with F aligned with one of the grid basis 
vectors, 3 of the 6 possible nearest neighbor directions 
will have a component parallel to F and so p = 1/2; for 
a square grid, p = 1/4. Thus, after one pass through the 
system, a randomly selected fraction p of the charges have 
advanced forward, while the rest remain in place. After 
a second pass through the system, a different randomly 
selected fraction p move forward. After many such passes 
one expects the system to be disordered. In fact, we find 
from simulations that, at T = 0, DDMMC disorders the 
ground state charge lattice for all F > F c . The random- 
ness of choosing moves, inherent in the DDMMC algo- 
rithm, is sufficient to disorder the moving system even as 
T ^ 0. 

We now consider the low T behavior of CTMC. For 
specificity we will consider the case of a triangular grid 
with a charge density of / = 1/25, and F = Fa\ aligned 
along one of the grid axes. We will study this particu- 
lar case in great detail in section ITTT1 Consider the limit 
T — ► 0, starting in the F = ground state and then 
increasing F, but with F < F c . The configuration of 
charges in the F = ground state is shown in Fig.|2K for 
a 25 x 25 grid. Since CTMC is a rejectionless algorithm, 
even when F < F c CTMC will make an excitation out of 
the ground state. However since AU > for this exci- 
tation, the time At for this excitation to occur diverges 
exponentially as T — ► 0. Conversely, once an excita- 
tion has been made, the very next move will be to relax 
the excitation back to the ground state, since this is the 
only move for which AU < 0; moreover, since AU < 0, 
this move takes place in an exponentially vanishing time. 
Thus alternating steps of CTMC will consist of displac- 
ing a randomly selected single charge and then moving it 
back. As T — + 0, the simulation clock time to be in the 
ground state grows exponentially large, while the clock 
time to be in the excited state gets exponentially small. 
The system therefore remains pinned in the ground state, 
with the time in the excited states contributing negligibly 
to any measured averages. 

Next, consider what happens when F > F c . Because 
of the rates of Eq. l|13f) . as T — > only moves with the 
smallest value of AUi a = AU m [ n can be selected; all other 
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FIG. 2: CTMC on a triangular grid with charge density / = 
1/25 at T — » and F > F c , with F parallel to the ai axis, 
(a) Ground state charge lattice for a 25 x 25 triangular grid. 
Numbers denote the locations of the charges in the ground 
state. The value of each number indicates the step on which 
that charge moves forward, (b) The change in interaction 
energy AH at each step as charges move forward. (•) are for 
a 25 x 25 grid and correspond to the moves in (a); (o) and 
(a) are the beginnings of similar sequences for 50 x 50 and 
100 x 100 grids. 



possible moves are exponentially suppressed. This results 
in the main difference between CTMC and DDMMC. In 
CTMC, as T — ► 0, motion is deterministic except for 
choosing randomly among moves with degenerate values 
of At/ m i n . Now consider starting in the F = ground 
state shown in Fig.^. All moves that advance a charge 
forward one grid spacing along di are equally likely with 
AC/min,! — F c — F < 0, while moves in all other directions 
are exponentially suppressed. Thus the first step will be 
to select any one of the N c charges at random and move 
it forward. On the second step, however, there are only 
two moves that have the new lowest At/mm, 2; these are to 
advance either the charge immediately in front of, or the 
charge immediately behind, the charge that moved in the 
first step. On the third step there are similarly only two 
moves with A{7 m i n ,3; advancing the charge immediately 
in front of, or immediately behind, the first two moved 
charges. The system proceeds in a similar manner until 
all charges in the same row parallel to F have moved for- 
ward. The next move will be to pick a charge at random 
in one of the two adjacent charge filled rows and move 



it forward, and then subsequently all charges in that row 
move forward one by one. In this manner, the rows of 
charges move one after another forward until the system 
has returned to the starting ground state, but with the 
entire charge lattice advanced by one grid spacing. In 
Fig-Et we label each of the ground state charges by the 
step number on which that charge moved forward in one 
particular pass of CTMC on a 25 x 25 size grid. The pat- 
tern described above is clearly evident. In Fig.[5p we plot 
the change in interaction energy, AH = AJ7 m i n ,„ +F, for 
each step n of this pass; note that AH is independent of 
the applied force F. The rough oscillation of AH with 
a period of n = Li/a$, with ciq = \f\ff the spacing be- 
tween the charges, reflects the row by row motion of the 
charges. 

Next we consider the timing of the above sequence of 
moves. From Fig. 03 we conclude that for each step n > 
1 of the above pass, t/ m i n ,n < fTmin.i- Hence the rate, 
Eq. I|13|) . to make any step n > 1 is exponentially larger 
than the rate to make the first step, n = 1. As T — > we 
conclude that the relative time spent in the intermediate 
states (i.e. the states after steps n = 1 . . . N c — 1) as 
compared to the time spent in the ordered ground state 
(prior to the first step and after step n = N c ) vanishes 
exponentially. According to the simulation clock time, all 
charges in the ground state charge lattice have advanced 
forward one grid spacing essentially simultaneously. This 
is the deterministic motion of the charge lattice that one 
would physically expect to find for F > F c as T — » 0. 

There is, however, one peculiar aspect to the above 
T — > dynamics. By the above arguments, the veloc- 
ity of the moving charge lattice will be proportional to 
the rate to make the initial first step. From Eqs. 1|13|) 
and (SJ this rate will be W to t = iV c Woe" AC/m,n < l/2T . 
Thus the T — ► velocity grows proportional to the num- 
ber of charges N c in the system. However this can be 
understood physically if one views motion on the dis- 
crete grid as being a representation for continuum mo- 
tion in a periodic potential. In this case, one should take 
W ~ e~ 2C/ °/ T , as in Eq. (J22J), where 2U is the maxi- 
mum to minimum energy difference of the potential. In 
the term W to t above, the factor Wo represents the rate for 
a particular charge to be excited out of the ground state, 
over the energy barrier 2Uq into the neighboring down 
stream potential minimum, thus lowering the energy of 
the system by AU min ^i. This rate becomes exponentially 
slow as T — > 0. Once this initial excitation has taken 
place, all the other charges follow, advancing forward in 
what may be regarded as an "avalanche". We call this 
an avalanche because all the other charges move forward 
in a time that becomes vanishingly small compared to 
the time to make the initial excitation. The factor N c in 
Wtot just reflects the N c possible sites at which the initial 
excitation that leads to the avalanche can occur. Thus, 
while Wtot grows as N c , it also vanishes exponentially as 
T — * 0, and so at T = the charges are always pinned, 
as is physically appropriate for F c < F < Uq. 

Several features of the expected behavior at finite T 
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can also be inferred from Fig.|5jD. Once a first charge in 
a given row has moved forward, the energy change for the 
other charges in the same row to move forward rapidly 
decreases. Consequently, once the first charge has moved 
forward, the remaining charges in that row rapidly follow 
forward. However, once a row has completely moved for- 
ward, the energy for the first charge in an adjacent row 
to move forward is not much lower than the energy for 
a first charge in any other row to move forward. Com- 
paring the values of AH. in Fig.^ for steps n = 1 and 
n = (Li/do) + 1, we estimate this energy difference as 
AE ~ 0.008. We therefore expect that once the temper- 
ature T become of the same order as this AE, coherence 
between moving rows will be lost. Avalanches will now 
consist of individual rows moving forward, but different 
rows will be uncorrelated. Consequently, the average ve- 
locity in this regime will scale proportionally to Li/oq 
(the number of sites in a given row for an initial excita- 
tion to occur) rather than N c . The details of this picture 
will depend on the specific correlations between charges 
within a given row, versus between rows, and this will be 
a subject of investigation in section ITTT1 



order parameter, 

$6 



,6iB, 



(27) 



In the above, the first sum is over all charges rii = 1, the 
second sum is over all charges rij = 1 that are nearest 
neighbors of rij, Zi is the number of such nearest neigh- 
bors, and 9ij is the angle of the bond from to rij with 
respect to the ai axis. Nearest neighbors are determined 
by Delaunay triangulationi^. 

We also measure several dynamical quantities. Let 



1 



c s<t 



(28) 



be the net displacement of the center of mass of the 
charges at time t of the simulation clock. The right hand 
side of the above is just the sum of charge displacements 
at each step s of the simulation that occurs before the 
simulation clock has reached time t, normalized by the 
total number of charges. The average velocity of the sys- 
tem is then just, 



C. Observables 

To determine the properties of our non-equilibrium 
steady states, we measure several static (time indepen- 
dent) and dynamic (time dependent) quantities. To de- 
termine structural properties, the main quantity of inter- 
est is the structure function, 



•S'(k) = -^r(n k n_ k ) 



where, 



n k 



(23) 



(24) 



is the Fourier transform of the charge distribution rij, 
and k is one of the allowed wavevectors in the first Bril- 
louin zone (shown in Fig.^D for the triangular grid). The 
corresponding real space correlation function is given by, 



C{m x ,m 2 ) = e- 2 ™ {rnikl+m2k2) 

ki,k2 



S(fci,*2) , (25) 



where we have expressed the positions and wavevectors 
k in terms of their coordinates and kn, as in Eqs. (J5J 
and J2J , in constructing the above Fourier transform. We 
will also consider the mixed correlation, 



C(fci,m 2 ) = ^2^ e 



S(k 1 ,k 2 ) . (26) 



The above quantities give information about the trans- 
lational order of the system. To investigate the orienta- 
tional order, we define the 6-fold orientational (hexatic) 



Rcmv ) Rcm(^eq) 



(29) 



where r is the total simulation clock time, and t cq is some 
initial time to allow for equilibration. In the analogy be- 
tween 2D charges and vortices in a superconducting film, 
the average charge velocity becomes the average voltage 
drop transverse to the direction of motion of the vortices. 

We also look at the fluctuations about the average cen- 
ter of mass position. If we define the fluctuation, after a 
time t, about the average center of mass position, 

SR(t] t ) = Rcm(f + t ) - Rem (to) - v avc i , (30) 

then we can define the diffusion tensor D(t) by, 

2B(t)t = N c (SK(t;t )5H(t;t ))t , (31) 

where the angular brackets in the above denote an aver- 
age over the parameter to during the course of a single 
simulation. In averaging over to , we restrict ourselves to 
non-overlapping intervals, i.e. to the values to — nt, for 
integer n, so as to reduce correlations among the different 
terms being averaged. If fluctuations about the center of 
mass are diffusive, then we expect D(t) to saturate to a 
constant as t increases. The factor N c on the right hand 
side of Eq. I|31|l ensures that D(i) approaches a size inde- 
pendent value in the liquid state, where the charges have 
only short ranged correlations. 

Although in our simulation we will use Eq. I|31|) to com- 
pute D(t), the diffusion tensor can also be expressed in 
its more familiar form, in terms of velocity correlations^. 
If we define the instantaneous fluctuation in velocity by, 



Sv(t) = v(t) - v a 



5H(At;t)/At , 



(32) 
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then 

lim D(i) = / dt(«5v(i)<Sv(0)) . (33) 
t^oo 2 J_ oc 

For a superconducting network, where vortex velocity is 
proportional to the voltage drop in the direction trans- 
verse to the vortex motion, D is a measure of the voltage 
fluctuations. 

In equilibrium, when F = v avo = 0, D/T is propor- 
tional to the charge mobility tensor by the fluctuation- 
dissipation theorem 1 ^. In the analogy to vortices in su- 
perconducting films, this is the linear resistivity of the 
film. In the driven state, with F = Fx, we will use the 
transverse component of the diffusion tensor, D yy , to test 
for the presence of transverse pinning. If D yy > 0, the 
center of mass is diffusing transversely to the direction 
of the average motion; application of a small transverse 
force SFy will cause the system to acquire a transverse 
component of the velocity, v y oc 5F. In analogy with 
equilibrium, we will assume that if D yy = 0, there will be 
no linear transverse response, i.e. v y /SF — > as 5F — > 0. 
This characterizes a transversely pinned stated. 

In CTMC, averages in the steady state are computed 
by the time integral in Eq. (|17|1 . However, the direct ap- 
plication of Eq. (|17|) would require the evaluation of the 
measured quantity after every single step of the simula- 
tion. For quantities involving lengthy calculation, such 
as S(k) and $6, this is not practical except for fairly 
short runs. Instead, we compute the time integrals for 
these quantities by a Monte Carlo integration-^, averag- 
ing them over N con fi g configurations sampled randomly, 
with a uniform distribution, over the simulation clock 
time interval (i cq , r), with i eq an initial equilibration time 
and r the total simulation time. In practice, we imple- 
ment this scheme as follows. We compute the average 
time interval between samplings, r' = (r — i oq ) / N con R g . 
Then, after a first sampling, we determine the time until 
the next sampling by drawing from an exponential dis- 
tribution with average time constant r'. This gives the 
correct sampling since, if t% < t 2 . . . < t n are the or- 
dered values of n independent and uniformly distributed 
random variables on a given interval, the probability dis- 
tribution for the distance tj+i —t% is exponential. We use 
typically N con a g ~ 10 3 to 10 4 in our simulations. 

III. RESULTS ON A TRIANGULAR GRID 

We now report our results for the case of charges on a 
triangular grid. The equilibrium, i.e. F = 0, behavior— of 
this system depends on the charge density /. For suffi- 
ciently dense / (but not too dense) , there is only a single 
first order melting transition at T m , from a pinned charge 
solid with long range translational order at low T, to an 
ordinary liquid at high T. For more dilute /, there are 
three phases: a low T pinned solid with long range trans- 
lational order, an intermediate T floating solid with al- 
gebraic translational order, and a high T liquid. In this 



work we will consider the charge density / = 1/25, which 
falls in the dense limit with a single first order equilib- 
rium melting temperature of T m ~ 0.0085. The dilute 
limit will be considered elsewhere. 



A. Phase Diagram 

We start by mapping out the T — F phase diagram for 
a 60 x 60 grid, with the applied force along the &i grid 
axis, F = Fx. We initialize the system in the ordered 
F = ground state, set F to the desired driving force, 
and then simulate the system for increasing values of T. 
By measuring the average interaction energy (H) , struc- 
tural properties such as S(k) and ($6), an d the average 
velocity v ave , we determine the phase diagrams shown in 
Fig-ED Our simulations consist typically of ~ 10 5 to 10 6 
passes through the system, depending on system size and 
parameters. 

Our results for the DDMMC dynamics are shown in 
Fig-Ot- As discussed earlier in section III B 31 at T = 
the system remains pinned (v ave = 0) in its equilibrium 
ground state for all F below a critical force F c = 0.0603; 
for all F > F c the moving state is a liquid^. For fixed 
F < F c , upon increasing T, the solid remains pinned with 
long range translational order, up until a value T P (F). It 
then enters a moving state with finite v ave . Over most 
of the phase diagram we find^ that this moving state is 
a liquid with a correlation length of order the spacing 
between charges do- Only in a very narrow region at 
low F and higher T do we find a structure that appears 
to be a moving smectic phase. We will discuss what we 
mean by a "smectic" phase in greater detail when we 
describe our results from CTMC. For DDMMC we have 
not investigated in any detail the stability of this small 
region of smectic phase with respect to increasing the 
system size, or with respect to cooling from the liquid. 
The lack of structure for almost all of the moving state, 
particularly at large F and small T, suggests that the 
DDMMC algorithm is indeed unphysical and unlikely to 
be a good model for continuum dynamics. We therefore 
focus on the CTMC algorithm for the remainder of this 
paper. 

In Fig.|3j3 we give the phase diagram for CTMC dy- 
namics. Again we find pinned, liquid, and smectic 
phases, but now the smectic persists over a wide region 
of the T — F plane, particularly at low T and large F. 
To illustrate the nature of order in each of these phases, 
we plot in Fig.0|the structure function 5(k) for several 
representative points in the phase diagram. In FigJ^Ji we 
see the sharp Bragg peaks with S(K) ~ N c , character- 
istic of the long range translational order in the pinned 
solid phase. The peaks are at the reciprocal lattice vec- 
tors of the ordered charge solid, and given by, 

K Pi,p 2 = kibi + k 2 b 2 

with *V = y> p M = 0,±l,±2 . (34) 
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FIG. 3: Phase diagram for a 60 x 60 triangular grid with 
charge density / = 1/25 as a function of temperature T and 
driving force F = Fx, (a) for DDMMC dynamics, and (b) for 
CTMC dynamics. "PS" stands for pinned solid. 



Fig. El 3 shows a roughly circular and finite peak, charac- 
teristic of short range translational order in the moving 
liquid phase. 

Figs.0Ji,d show the moving smectic at small and large 
driving forces, respectively. Consider first the large drive 
case in Fig.^Ji. The peaks along the ki axis [k\ — 0) at 
(fcijAfe) = (0,±l/5) and (0,±2/5) (see Fig.Ht for the k- 
space geometry) are sharp Bragg peaks with S(K) ~ N c . 
This indicates that if one averages the particle density in 
the cii direction (fci = 0), the resulting density is peri- 
odic in the &2 direction with a period of 5 grid spacings; 
the particles arc therefore moving in periodically spaced 
channels oriented parallel to the driving force F = Fa^. 
Next, note that the peaks at finite k\ = ±1/5, ±2/5 ap- 
pear to be sharp, i.e. only one grid spacing wide, along 
the k\ direction. Such <5-function like peaks in the k\ 
direction indicate that the particles are periodically or- 
dered within each smectic channel with a period of 5 grid 
spacings. The finite width of these peaks in the k-i di- 
rection indicates that the ordered smectic channels are 
randomly displaced with respect to each other, with a 
finite correlation length £j_ proportional to the inverse 
width of the peak. Comparing Fig.^Ji with Fig.^Ji, we 
see similar features at the smaller drive F, only the peaks 
at finite k\ are now sharper, with a narrower width in 
the ki direction. In the next two sections we will con- 
sider these features of the smectic phase in much greater 
detail, studying the scaling behavior and stability of the 
smectic as the system size increases. 

Finally we consider the nature of the melting transition 
T m (F) from the smectic to the liquid, and the unpinning 
transition T p (F) from the pinned solid to the smectic. 
In Fig.[SjL— d we plot the average interaction energy per 
site, E = (H/N) vs. the number of simulation passes 
through the grid. Each point represents an average over 
3200 successive passes. Fig.|SJi shows results at F = 0.02, 
T = 0.0042, just at the melting transition T m (F). We see 
that the energy takes a discontinuous jump as the sys- 
tem makes the transition from smectic to liquid. Melting 
of the smectic is therefore like a first order phase tran- 
sition. Fig.Eb, shows results at F = 0.02, T = 0.0022, 



just above the unpinning transition T p (F). We see that 
the energy fluctuations form a set of plateaus, with a 
very long period of fluctuation. The lowest plateau cor- 
responds to the ordered F — ground state with a value 
Eq = 0.03495736. The higher energy plateau corresponds 
to having some fraction of adjacent smectic channels (i.e. 
charge filled rows) advanced one grid spacing parallel to 
F, so that the system looks locally like the ground state, 
but with one pair of domain walls parallel to F. 

As T increases above T p , Figs.[3;,d show that the rate 
of fluctuations increase, and plateaus of additional energy 
values appear. The higher energy plateaus correspond to 
having more than one pair of domain walls in the other- 
wise ordered system. This behavior may be understood 
by considering the results shown in Fig.[2x For a driving 
force of F = 0.02, the thermal energy needed to excite 
a pair of adjacent particles in a given row to move for- 
ward is AU = AJ7 m i„4 + Ai7min,2 » 0.048; however the 
energy to move each remaining particle is A£/ m i n ,n < 0. 
Thus, at low T, the excitation of an initial pair forward 
will trigger the remaining particles in that row to move 
forward almost instantaneously. The rate for the initial 
pair excitation goes as, W ~ e -A(7/2T^ van ishing expo- 
nentially as T decreases. The rate of energy fluctuations 
decreases accordingly. At low T, once a given row has 
moved forward, the next most favorable excitation is to 
move an adjacent row of charges forward (see discussion 
at the end of section III B 3|) . The system thus consists 
of a single pair of domain walls in the otherwise ordered 
ground state; the distance between the domain walls in- 
creases as more adjacent rows move forward. Such states 
give the higher of the two energy plateaus in Fig.JSjD. 
As T increases, there becomes a non-negligible probabil- 
ity to have a pair excitation in a non-adjacent row of 
charges. Now the system can develop more than a sin- 
gle pair of domain walls, leading to the additional high 
energy plateaus of Figs.|3J:,d. 

The above arguments suggest that T P (F) may not be 
a true phase transition. Since the above rate W is finite 
at any T, but grows vanishingly small as T — > 0, the ob- 
served T p (F) may just result from 1/W growing larger 
than the longest simulation time we can carry out. As F 
decreases, it will becomes necessary to excite three, then 
four, then more, particles forward in a given row, before 
one reaches the condition that At/ m i n ,n < triggering the 
remaining particles in the row to move immediately for- 
ward (see Fig.|5]D). Thus we expect that W will decrease, 
and the observed T P (F) will increase, as F — > 0. Note 
that the graphs in Fig.[S^— d are plotted versus the num- 
ber of simulation passes rather than the simulation clock 
time. They therefore reflect the amount of actual compu- 
tation needed to observe fluctuations of the system. The 
decrease in fluctuation rate observed as T decreases for 
fixed F = 0.02 (compare Fig.[3i to Fig.JSJj) results from 
the decrease in probability to move the second particle 
of an excitation pair forward, before the first particle has 
had a chance to fall back into place, rather than being 
directly due to the overall exponential decrease with T 




FIG. 4: S(k) for several points in the CTMC phase diagram of Fig.|3Jj- (a) F = 0.02, T = 0.0014 in the pinned solid; (b) F0.02, 
T = 0.008 in the moving liquid; (c) F = 0.02, T = 0.003 in the moving smectic; (d) F = 0.10, T = 0.004 in the moving smectic 
at higher drive. The bottom row shows intensity plots of the corresponding graphs in the top row. The peak S(k = 0) = N c is 
removed to give greater contrast to the other peaks. 



of all single particle rates as in Eq. 1|13|) • 

Finally we note that, for finite system size, the moving 
smectic state is the true stable steady state of the sys- 
tem. In Fig.03 we plot the average interaction energy per 
site E versus simulation clock time t, for the parameters 
F = 0.03, T = 0.003, which lies just immediately below 
the melting transition T m (F) (see Fig.^). Comparing 
with Fig.[5Ji, we see clearly that the system has occa- 
sional fluctuations into the liquid phase, as indicated by 
the large brief spikes in energy. Such fluctuations are ex- 
pected for a finite size system near a first order phase 
transition. The fact that the system returns to the smec- 
tic state, after such a liquid fluctuation, indicates that 
the smectic is indeed the stable steady state. We have 
also succeeded in cooling into the smectic state from the 
disordered liquid, and in entering the smectic from the 
liquid by increasing F at temperatures above the mini- 
mum of the T m (F) transition boundary. 



B. Smectic Phase - High Drive 

In the next two sections we explore in detail the nature 
of ordering in the moving smectic state, and its stability 
as the system size increases. We start here by considering 
the smectic in the high drive case at F = 0.1, T = 0.004, 
corresponding to Fig.^Ji. If the system has true long 
range smectic order, we expect the peaks in S(k) along 
the ky, axis {k\ = 0) to be true Bragg peaks, with a 
height that scales as the system area. In Fig.|Hlwe plot 



the height of the peak S(Koi), versus system length L, 
for systems of size Lx L. The straight line on the log- log 
plot has a slope s ~ 1.99 giving good agreement with the 
~ L 2 behavior expected for long range smectic order. 

Next we consider the ordering within the smectic chan- 
nels. If charges have long range order within each indi- 
vidual channel, we expect the peaks in S(k) at fei = 
±1/5, ±2/5 to be <5-function like in the k\ direction. If 
the channels have only short range correlations between 
them, the width of these peaks will remain finite in the 
&2 direction. We therefore expect that the heights of 
these peaks at finite k\ should scale as the system length 
L\ in the a\ direction. In Fig. [5] we plot the height of 
two of these peaks, 5(Kn) and S , (K 2 i) (see Eq. (|3*4"Jl for 
our notation labeling the reciprocal lattice vectors K) 
versus system length L, for systems of size L x L. The 
straight lines have slopes s — 1.15 and 0.96 respectively, 
in reasonable agreement with the ~ L behavior described 
above. 

To further illustrate the above results, we plot in Fig.0 
profiles of S(k) along different paths through the first 
Brillouin zone, for various L x L system sizes. Fig.EJi 
shows S(k) versus fci for fixed k% = 1/5. The loga- 
rithmic vertical scale, varying over five orders of mag- 
nitude, indicates how sharply the peaks are confined to 
the values k\ — 1/5,2/5; however S(k) appears to de- 
crease continuously as one moves away from the peak 
values. Fig.03 shows the peaks indicating the smectic 
order, i.e. S(\t)/ fL 2 versus fc 2 for fixed fci = 0. We see 
that the peaks, scaled by N c = fL 2 , all have the same 
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FIG. 5: (a)-(d): Average interaction energy per site E = 
(Tt/N) vs. the number of simulation passes through the grid. 
Each point represents an average over 3200 successive passes; 
(a) F = 0.02, T = 0.0042, just at the melting transition 
T m (F); (b) F = 0.02, T = 0.0022, just above the unpinning 
transition T P (F); (c) F = 0.02, T = 0.0026 and (d) F = 0.02, 
T = 0.0032, moving away from the unpinning transition, (e) 
E vs. simulation clock time t, for F = 0.03, T = 0.003, 
slightly below the melting transition. 



1000 



100 



00 



10 









T = 


0.004 




- F = 


0.1 




f = 


1/25 








^ii m 















20 



40 



60 80 100 
L 



FIG. 6: Scaling of peak heights S(K) vs. system size L for 
the smectic phase at F = 0.10, T = 0.004. Straight lines 
indicate good power law fits, S(K.) ~ L s , with s ~ 1.99 for 
Koi, s ~ 1.15 for Kn, and s ~ 0.96 for K21. 



height for the different L, in agreement with the scaling 
seen in Fig. [5] The logarithmic vertical scale, dropping 
five orders of magnitude as one moves a single grid point 
in fc-space away from the peaks at ki = 1/5,2/5, shows 
that these are indeed sharp Bragg peaks. In Figs.[7|:,d 
we show the peaks at finite ki, plotting S(k)ao/L versus 



fc 2 for fixed ki — 1/5 and k% = 2/5 respectively. We 
see that these profiles, when scaled by 1/L, collapse rea- 
sonably well to a common curve for the different sizes 
L, for all values of ki. This is in agreement with the 
scaling found in Fig.|H| and suggests that S(k) is indeed 
(5-function like in k\ = 1/5, 2/5 for all ki- 

The finite widths in the ki direction of the peaks in 
5(k) at k\ — 1/5,2/5, which do not narrow as L in- 
creases, indicate that the ordered smectic channels have 
only short range correlations between them. To see this 
explicitly, consider the Fourier transform of the charge 
density in each row of the grid at the wavevector cor- 
responding to the periodic ordering within the smectic 
channels, i.e. k = (l/5)bi, and compute the correlations 
of this Fourier amplitude between different rows. This is 
given by the correlation function C(fci = 1/5, mi), de- 
fined in Eq. JUJ. In Fig-EK we plot C(h = 1/5, m 2 ) 
versus 1112 for a 60 x 60 system. The correlation has 
peaks at values mi = 5n, n integer, and is essentially 
zero in between, indicating that the particles flow in pe- 
riodically spaced channels, and that the channels have a 
width of a single grid spacing. The exponential decay of 
the peak heights indicates the short range correlation be- 
tween particles in different smectic channels. In Fig.[SJi> 
we plot only the peaks of C(fci = 1/5, mi), but for dif- 
ferent system sizes L x L. The curves for different L lie 
almost on top of each other and decay to zero, indicating 
a finite, size independent, correlation length £j_ trans- 
verse to the direction of the applied force F. To estimate 
£_l we fit to a simple periodic exponential, 

C(ki = 1/5, ma) ^ A (V™ 2 ^ + e -( L - m *)/^ 

(35) 

and get values in the range £j_ ~ 7.0 ± 0.5 as L varies 
from 60 to 100. Thus correlations extend only slightly 
beyond nearest neighbor channels (which are separated 
by 5 grid spacings). 

Next we compute the correlations within individual 
smectic channels. In Fig.|5Ji we plot the real space corre- 
lation parallel to the driving force, C(mi, m 2 = 0) versus 
mi, for a system of size 60 x 60. Again we see sharply de- 
fined peaks at mi = 5n, n integer, corresponding to the 
periodic spacing of particles within the channel. More- 
over the height of these peaks decays only slightly to a 
large finite value as mi — > L/2, as one would expect for 
long range order. However, when we plot in Fig.|5jD the 
height of these peaks for different values of L for system 
sizes L x L, we now see behavior inconsistent with long 
range order. The value of C(mi,0) at any given value 
of mi decreases as L increases; the magnitude of this 
decrease from unity is proportional to L. Rather than 
indicating long range order, such behavior is consistent 
with a very dilute but finite density of order destroying 
defects; when L is small compared to the average spacing 
between defects, then the probability to have a defect in 
the system will be proportional to L, resulting in a de- 
crease in the correlation proportional to L. 

Indeed, since the smectic channels are essentially de- 
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FIG. 7: Profiles of S(k) in various directions, for different 
system sizes L, for the smectic phase at F = 0.10, T = 0.004. 
(a) S(k) vs. ki for fixed k 2 = 1/5; (b) S(k)/fL 2 vs. fc 2 for 



fixed ki = 0; (c) S(\t)ao/L vs. 
S(lt)ao/L vs. &2 for fixed fci = 
scale in (a) and (b). 



fca for fixed fci = 1/5; (d) 
2/5. Note the logarithmic 
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FIG. 8: Transverse correlation function C(ki = 1/5, m.2) vs. 
m 2 at F = 0.1, T = 0.004. (a) C(ki = 1/5, m 2 ) for all in- 
teger values 7Ti2 for the single size L = 60; the dashed line 
highlights the decaying envelop of the peaks, while the solid 
fine interpolates between the data points, (b) Peak values 
of C(fci = 1/5,7712) at values 7712 = 5n, n integer, for differ- 
ent sizes L; solid lines are fits to the periodic exponential of 

Ea. 1551 . 



coupled from each other (as illustrated in Fig.|SJ|, each 
channel can be thought of as an independent one dimen- 
sional system. Although the charges in the channel have 
a bare long range logarithmic interaction, the uncorre- 
cted motion of charges in neighboring smectic channels 
will screen this log interaction, converting it to an ef- 
fective interaction that is finite ranged. In equilibrium, 
such a one dimensional system must be disordered at any 
finite T, and we expect that the same will be true of a 
driven steady state. To test this we perform independent 
CTMC simulations of a one dimensional (ID) lattice gas 
of particles with average spacing 5 and a nearest neighbor 
harmonic interaction with a spring constant k. Carrying 
out simulations for the same system sizes L as in Fig.^, 
we adjust k to get the best fit to the correlations found in 
the original two dimensional system. This gives a reason- 
able value of k — 0.0505, and our ID results are shown as 
the open symbols and dashed lines in Fig.|5jD. We see that 
the agreement is very good; the small deviations that ex- 
ist are presumably due to the small but finite coupling 
between neighboring smectic channels that exists in the 
original 2D model. Having found k, we can now simulate 
the ID model for much larger L, to see the exponential 
decay of correlations in the model and to determine the 
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FIG. 9: Longitudinal correlation function C(mi,m,2 = 0) 
vs. mi at F = 0.1, T = 0.004. (a) C(rai,0) for all inte- 
ger values mi for the single size L — 60. (b) Peak values of 
C(mi,0) at values mi = 5n, n integer, for different sizes L 
(solid symbols) ; solid lines are fits to the periodic exponential 
of Eq. 11361 . Open symbols and dashed lines are fits to a one 
dimensional model (see text). 



correlation length We find for the ID model, £|| ~ 86. 
For comparison, we can fit the data from the original 2D 
model to a periodic exponential, 

C(mx,m 2 = 0) ~ ci +c 2 (e^™ 1 ^" + e - (L - mi)/ ^) , 

(36) 

where c\ — 1/5 is the average density of charges in 
a smectic channel, and C2 = 4/5 is chosen so that 
C(0,0) = 1. The resulting fits are shown as the solid 
lines in Fig.|U>, and give the values £|| ~ 83,80,78 for 
sizes L = 60, 80, 100, in good agreement with the ID 
model. We conclude that the smectic phase at high 
drive consists of weakly coupled channels, characterized 
by a small transverse correlation length £j_ . Within each 
channel particles have only finite range correlations, but, 
for the case considered above, this longitudinal corre- 
lation length is comparable to the size of the system, 
£|| ^ L 3> so that the particles in a given channel 
appear to be ordered. 

We can next ask what happens if the system length 
parallel to the driving force F increases to be larger than 
the longitudinal correlation length, L\ > £«. Increasing 
the system to size L\ x L2 = 500 x 60, so that this condi- 
tion is met, we found in Ref. that the smectic phase at 
F = 0.1, T = 0.004, becomes unstable to the liquid; for 
this system with bigger L\, the smectic remains stable 
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FIG. 10: Correlation function (7(mi,77l2 = 0) vs. mi at F = 
0.1 and various T, for a 60 x 60 system. Only the peak values 
at mi = 5n, n integer, are shown. Solid lines are fits to an 
appropriate periodic exponential, as in Eq. 11361 , excluding the 
initial point at mi = from the fit. 



only at lower T such that the condition L\ < £11 is again 
obeyed. To illustrate this point further, we carry out 
simulations for a system of size 60 x 60 at F — 0.1, but 
increasing T so as to cross the melting line shown in the 
phase diagram of Fig.|3|3. In Fig.^H we plot the result- 
ing correlation functions C(mi,m2 = 0) versus mi for 
several different temperatures. We see clearly the tran- 
sition from smectic to liquid at a temperature between 
0.005 and 0.006. To get an estimate of the longitudi- 
nal correlation length £|| we fit the data of Fig.^5] to 
a periodic exponential, as in Eq. (|36|) . For the smectic, 
T < 0.005, we set ci = 1/5 as appropriate for the av- 
erage density of charges in a smectic channel. For the 
liquid, T > 0.006, we set c\ — 1/25, appropriate for the 
average density of charges in the liquid. In both cases we 
find better results when we exclude the initial point at 
mi = 0, C(0, 0) = 1, from the fit; we therefore keep C2 as 
a free fitting parameter. The resulting fits are shown as 
the solid lines in Fig.^H The values of £|| obtained from 
these fits are then plotted versus 1/T in Fig.lTTk The 
dashed straight line on the plot indicates an Ahrenius 
form, £|| ~ e T °/ T , for the divergence of £|| in the smectic 
as T — > 0. We see that melting occurs when £|| ~ 30, i.e. 
roughly half the system length. We conclude that the 
smectic is only stable when £|| > L/2. When the corre- 
lation length becomes smaller than this, and one would 
expect particles within the smectic channels to disorder, 
the entire smectic structure becomes unstable to a liquid. 
Since £|| diverges rapidly as T — > 0, however, one should 
expect to see the smectic phase in any finite size system, 
at sufficiently low temperature. 

Finally, we can also estimate the transverse correlation 
length £x< F° r the smectic we use an analysis of C(ki = 
l/5,m2), similar to that of Fig.|HjD, to determine £x- F° r 
the liquid, since there is no periodic ordering, we use 
an analysis similar to that of Fig.^| applied to the real 
space correlation C{x — 0,y) = C(m\ = —1712/2,9712) 
versus y = (v / 3/2)m 2 . Note that since the direction a 2 
is not orthogonal to hi (see Fig-QJi) it is necessary to 
use the argument mi = — 7712/ 2 in order to measure a 
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FIG. 11: (a) Longitudinal correlation length £m vs. 1/T, and 
(b) transverse correlation length £± vs. T, at F = 0.1 for a 
60 x 60 system. 



strictly transverse correlation. Our results for £j_ versus 
T are shown in Fig.lllb. Unlike the rapid rise in £|| as T 
decreases, we see only a small increase in £j_ as T — ► 0. 
In the liquid, £j_ and £|| are comparable. 



C. Smectic Phase - Low Drive 

We now consider the structure of the smectic in the 
limit of smaller driving forces, in particular the case F = 
0.02, T = 0.003, shown in Fig.BJ;. In FigJH we plot 
profiles of S(k) along different paths in the first Brilloun 
zone, for different system sizes L x L. Comparing with 
the analogous Fig .0 for the high drive case, F = 0.1, we 
see that the peaks are now more sharply defined. The 
peaks along the k 2 axis at k\ = in Fig. ll2b continue 
to look like Bragg peaks, being only one grid point wide 
and with heights scaling as L 2 , thus indicating long range 
ordering into smectic channels. However the peak heights 
at finite ki = 1/5, 2/5 in Figs.ll2b.d no longer appear to 
scale ~ L as do the corresponding peaks in Fig.0 

In Fig.ll3b we plot the heights of various peaks S(K) 
versus L for various system sizes L x L. While the 
smectic peak S'(Koi) scales ~ L 2 as expected, we find 
S(K-u) ~ L 13 , more divergent than the ~ L behav- 
ior found at higher drive. This suggests the possibility of 
longer range, perhaps algebraic, correlations between the 
different smectic channels. For algebraically diverging 
peaks, however, we expect that not only the peak height 
must scale, but the entire peak profile should scale. The 
expected scaling relation is^X, 



S(k 1 ^l/5,k 2 )^L 1 - 3 f(k 2 L) 



(37) 



where f(x) is a scaling function. In Fig.ll3b we test this 
scaling prediction by plotting S(k\ = 1/5, k 2 )/L 13 ver- 
sus k 2 L for different sizes L. We clearly do not find the 
scaling collapse expected for algebraic correlations. 

To explain the above behavior, we consider the trans- 
verse correlation function C{k\ — 1/5, m 2 ), which we plot 
versus m 2 for different system sizes L x L in Fig. ll4h ,. 




FIG. 12: Profiles of S(k) in various directions, for different 
system sizes L, for the smectic phase at F = 0.02, T = 0.003. 
(a) S(k) vs. fci for fixed k 2 = 1/5; (b) S(k)/fL 2 vs. k 2 for 
fixed ki — 0; (c) 5 , (k)ao/I/ vs. k 2 for fixed fci = 1/5; (d) 
S(\t)ao/L vs. k 2 for fixed ki — 2/5. Note the logarithmic 
scale in (a) and (b). 
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FIG. 13: Smectic phase at low drive, F = 0.02, T = 0.003. 
(a) Scaling of peak heights S(K) vs. system size L. Straight 
lines indicate good power law fits, S'(K) ~ L s , with s ~ 2.0 
for Koi, s ~ 1.3 for Kn, and s ~ 0.87 for K21. (b) Attempted 
scaling collapse of S(ki = l/5,fc2)/I/ vs. k^L. 



The solid lines are fits to the periodic exponential of 
Eq. 1(231, and give the values £_L ~ 29.6,30.8,29.7 for 
sizes L = 80, 100, 140, respectively. Our results thus con- 
sistently indicate short range order between the smec- 
tic channels, with a finite transverse correlation length 
~ 30. The absence of the expected ~ L scaling in 
Figs.ll2b.d is then a finite size effect due to the corre- 
lation length £j_ being comparable to the system length 
L. 

We similarly estimate the longitudinal correlation 
length by plotting C(mi, iri2 — 0) versus mi, for different 
system sizes L x L, in Fig.ll4b. Fitting to the periodic 
exponential of Eg. 1|36[1 . with c\ — 1/5 and C2 = 4/5, 
we find £|| ~ 4700. We conclude that the smectic phase 
at small drive is qualitatively the same as that at large 
drive, except for having larger, but still finite, correlation 
lengths. 

Finally we consider behavior as the driving force F 
varies. In Fig. ll5b we plot the transverse correlation 
function C(ki = 1/5, 7712) versus mz, for T = 0.003 
in a 60 x 60 system, for various values of F from the 
low drive case considered above, F — 0.02, to the high 
drive case considered previously, F = 0.1. Solid lines are 
fits to the periodic exponential of Eq. i|35|) . In Fig.ll5b 
we plot the corresponding longitudinal correlation func- 
tion C(mi,TO2 = 0) versus mi. Solid lines are fits to 
the periodic exponential of Eq. (|36|l . with c\ = 1/5 and 
C2 = 4/5. From these fits we estimate the transverse 
and longitudinal correlation lengths, £j_ and £||, which 
are plotted in Figs.ll6b,.b. We see that the transverse 
correlation length £j_ grows as F decreases below 0.05. 
For F > 0.05, £j_ levels off to a constant, £j_ ~ 7. In fact, 
for F > 0.05, the transverse correlation C(k\ — 1/5, 1712) 
shown in Fig.ll5b is completely independent of F. The 
longitudinal correlation length £|| grows exponentially 
as F decreases below 0.04, has a shallow minimum at 
F = 0.05 (presumably related to the minimum in the 
melting line T m (F) that occurs nearby) and then satu- 
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FIG. 14: Smectic phase at low drive, F = 0.02, T = 0.003, for 
various system sizes L x L. (a) Transverse correlation func- 
tion C(fci = 1/5, mj) vs. 77i2, and (b) longitudinal correlation 
function C(mi,m,2 ~ 0) vs. mi. Solid lines are fits to peri- 
odic exponentials as in Eqs. I|35^ and l|36|l and determine the 
correlation lengths £± and £y . 



rates to a constant above F = 0.06. The longitudinal 
correlation C(n%i,m2 = 0) shown in Fig. 115b is indepen- 
dent of F for F > 0.06. These results strongly suggest 
that no new behavior will be seen by increasing F to 
even larger values, and we have verified this explicitly by 
simulating up to F — 2.0. 

The above cross over point between low and high 
drives, and its value F cv ~ 0.05, can be understood from 
Fig-Eb- Consider the system in its F = ground state. 
The interaction energy to move a given charge forward 
parallel to F is AHi ~ 0.0627; the rate to make this move 
is W = Woe~^ A ' Hl ~ F ^ 2T '. Once this charge has moved 
forward, the interaction energy to move the neighboring 
charge in the same row forward is AH.2 — 0.0270; the rate 
to make this move is Wf = Woe~( An2 ~ F ^ 2T . This needs 
to be compared against the rate for the first charge to 
move back to its original position, W\> = Woe^ Ani ^ F ^ 2T . 
The ratio of these last two rates is, 



3 (2F-A"Hi-AH 2 )/2T 



(38) 



and so the two rates are equal when F cr = (A7ii — 
A7i2)/2 = 0.045, which agrees with our observations in 
Figs. E>1 and [Tj3 When F > F cr , the probability that the 
neighboring charge moves forward along with the first 
charge is larger than the probability that the first charge 
falls back into place. Once the second charge moves for- 
ward, the remaining charges in the row will follow suite. 
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FIG. 15: Smectic phase at T = 0.003, for various driving 
forces F in a 60 x 60 size system, (a) Transverse correlation 
function C(k\ = 1/5, 7712) vs. 7712, and (b) longitudinal corre- 
lation function C(nii,m2 = 0) vs. mi. Solid lines are fits to 
periodic exponentials as in Eqs. 1351 and I36H and determine 
the correlation lengths £± and f y . 
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FIG. 16: (a) Transverse correlation length £y and (b) longitu- 
dinal correlation length £± vs. F, at T = 0.003 for a 60 x 60 
system. 



Therefore when F ^S> F CT , virtually all moves are those 
which advance a charge in the direction of the applied 
force F; charges in the smectic channels move continu- 
ously forward row by row. For F < F CT , charges which 
advance forward parallel to F will more often than not 
fall backwards to their original position on the next move. 
The system spends finite time with no net motion, in be- 
tween randomly occurring avalanches that advance an 
entire row of charges forward. The result is a stick-slip 
type of motion. 

The above scenario is illustrated in Fig. 1171 where we 
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FIG. 17: Center of mass displacement parallel to the driving 
force for T = 0.0022, L = 60. N c X cnl vs. simulation clock 
time t for (a) F = 0.02 and (b) F = 0.05. Insets show an 
expanded picture on a short time scale. 



plot the center of mass displacement parallel to F times 
the number of charges, N c X cm , versus the simulation 
clock tme t. We show results for T — 0.0022, slightly 
lower than the temperature 0.003 considered above, for 
a system of size 60 x 60. Fig.[T3i, for F = 0.02 < F cr , 
shows the step like advancement forward of the system, 
characteristic of stick-slip motion. The inset shows an 
expanded scale for short time. The height of each step 
is exactly 12 grid spacings, corresponding to the advance 
of all 12 charges in a given smectic channel. Along the 
plateau of each step we see motion one grid space for- 
ward, followed by one grid space backwards, with no net 
motion. In Fig.ll7b we show results for F = 0.05 > F CT . 
We see that motion is perfectly linear in time. The inset 
shows that the system moves smoothly forward, with one 
row advancing immediately after another. 



D. Liquid Phase 

We now briefly consider the liquid phase. The liquid 
phase shown in Fig.^J), at F — 0.02, T = 0.008, appears 
fairly structureless. However as T decreases, and the fi- 
nite correlation lengths grow, local structure develops. 
As discussed in section UlI Bl the melting line T m (F) de- 
creases as the size of the system increases. Although, for 
L = 60, T m (F) always lies above T = 0.003 (see Fig.|3>), 
when L increases, the minimum of T m (F) near F ~ 0.03 
dips below 0.003. In Fig. ^| we plot the structure func- 
tion 5(k) for a system of size 120 x 120 at T = 0.003 and 
driving force F = 0.05. Although one sees prominent 
peaks at the reciprocal lattice vectors corresponding to 
the F = ground state, the system is in a liquid state 
with short range translational order. The heights of the 
peaks are small compared to those in a more ordered 
(i.e. smectic or solid) state, and as L increases, the peak 
heights stay finite rather than diverging with system size. 

Next we consider how behavior in the liquid varies with 
the driving force F. In Fig.lliJb we plot the longitudinal 
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FIG. 18: Structure function S(k) for system of size 120 x 120 
in the liquid state at T = 0.003 and F = 0.05. 
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FIG. 19: (a) Correlation lengths £± and £m vs. T for a system 
of size 120 x 120 in the liquid state at T = 0.003. (b) Orien- 
tational order parameter | ( *3?e ) | vs. F for several systems of 
size L x L at T — 0.003. The larger valued data points are 
for the smectic phase; the lower data points are for the liquid 
phase. 



and transverse correlation lengths £j_ and £ii versus F. 
We obtain £|| and £j_ by fitting to the correlation func- 
tions C(mi, TO2 = 0) and C(x — 0, y) in the same way as 
we have done earlier in constructing Fig.^] We see that 
£j_ and ^ increase with increasing T, and that £,± ~ 2£|| . 
That order extends further in the transverse than the lon- 
gitudinal direction can be seen by noting that the trans- 
verse peaks in 5(k), shown in Fig. 1181 are larger than the 
peaks with a longitudinal component. The same obser- 
vation was made in our earlier work 8 for a much larger 
system at higher driving force. 

We also investigate orientational order in the liqud. In 
Fig. 119b we plot the absolute value of the 6-fold orienta- 
tional order parameter |($6)| versus F at T = 0.003, for 
several different system sizes L x L. Depending on the 
value of L and the corresponding value of T m (T; L), the 
system is either in the smectic state (with a high value of 
| ($6) I) or m the liquid state (with a low value of |($e)|). 
Even in the liquid | ($6) I is finite because of the 6-fold 
rotational symmetry of the underlying triangular grid. 
We see that, as F increases in the liquid, |($e)|> like £i 
and £|| , increases. Thus both translational and orienta- 
tional order in the liquid increase as the driving force F 
increases. 



FIG. 20: Average velocity ii ave x parallel to the driving force 
F for a 60 x 60 system, (a) v^ vcx vs. 1/T for fixed F — 0.1 
in the high drive limit, (b) v avcx vs. F for fixed T = 0.003 
in the smectic. The dashed lines indicate the exponential 
dependence of v mcx on F and 1/T. 



E. Dynamics 

The preceding sections have dealt with the structural 
behavior of the driven system. In this section we con- 
sider some of the dynamical behavior. In Fig.[5n| we 
plot results for the average velocity n av „ parallel to the 
driving force F = Fx, for a 60 x 60 size system. In 
Fig.l2Hi we show results for fixed F — 0.10 > F cr , in 
the high drive limit, versus 1/T. The dashed line for 
1/T m < 1/T shows the exponential dependence on 1/T 
in the smectic phase, u avel ~ e T °^ T . Fitting to this form 
we find To ~ 0.02. We can understand this value as 
follows. As discussed in section IIII CI at such high F 
virtually all moves in CTMC result in the advance of 
a charge forward; the charges in the smectic channels 
move steadily forward one channel at a time. The aver- 
age velocity is set by the rate for the charges in a given 
channel to move forward, which in turn is set by the rate 
for the first charge in the channel to move forward (all 
other charges in that channel moving forward on a much 
more rapid time scale). The rate for the first charge 
in a channel to move forward is ~ e -A[/ mi „i/2T^ wnere 
-A[/ min i = F - AHi, with AHi ~ 0.06 from Fig.03. 
Thus we have T = (F - AHi)/2 = 0.02. 

In Fig.l20b we show results for v avcx vs. F, for fixed 
T = 0.003 < X" m , in the smectic phase. The dashed 
line shows the exponential dependence on F in the high 
drive limit, F > T cr ~ 0.045, v avcx ~ e F/F °. From the 
preceding discussion, we expect To = 2T = 0.006, and 
we find this value gives an excellent fit to the data. 

We have also considered the dependence of the average 
velocity on the system size. In Tabled we list the values 
of Wave a: for various system sizes, with different aspect 
ratios Li/T 2 , at F = 0.10, T = 0.004, in the high drive 
smectic. In agreement the discussion at the end of section 
III B 31 we see that u avel scales roughly proportional to the 
length of the system L\ parallel to the driving force. 

Next we consider the diffusion of the center of mass 
about its average motion. To compute D(t) we need 
to compute the correlation between states of the sys- 
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TABLE I: Average velocity ti avcl for various system sizes on 
a triangular grid at F = 0.10, T = 0.004, in the high drive 
smectic. 
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FIG. 21: Center of mass diffusion constants (a) D yy and (b) 
D xx vs. time t for the smectic phase in the low drive limit, 
F = 0.02, T = 0.002, for a system of size 60 x 60. That 



indicates the system is transversely pinned. 



tem separated by time t. Since we are interested in 
the long t limit, computing D(<) accurately thus requires 
much longer simulations than were needed to compute 
the structural (equal time) correlations. We therefore 
present results only for several typical cases. In Fig.EU 
we plot D yy (t) and D xx (t), defined by Eq. J^SJ, versus 
the simulation clock time t, for the smectic in the low 
drive limit of F = 0.02, T = 0.002, for a 60 x 60 size sys- 
tem. We see that D yy decays to zero, indicating that the 
system is transversely pinned. D xx saturates to a finite 
constant as t increases, indicating a random walk motion 
about the average center of mass position. In Fig.1221 we 
similarly plot D yy and D xx for the smectic in the high 
drive limit of F = 0.05, T = 0.0022. Again we see that 
D yy — > and the system is transversely pinned, while 
D xx saturates to a finite constant. In Fig. 1231 we plot 
D yy and D xx for the liquid at F = 0.10, T = 0.006, just 
above the melting transition. In this case both D yy and 
D xx approach finite constants as t increases; as expected, 
the liquid is not transversely pinned. 

Since both D xx and n avcl approach constants in the 
long time limit, a convenient measure of the strength 
of fluctuations about the average motion is given by 
D xx /v avcx . For the low drive smectic of Fig. 1211 we find 
D xx /v avcx ~ 40. This is consistent with our inter- 
pretation of this region as being one of stick-slip mo- 
tion. In this case we expect that the motion of rows of 
charges forward will constitute a Poisson process with 
avalanches occurring at a rate A. At each avalanche 
n r y/fL charges move forward, where \fjL is the number 
of charges in a given smectic channel, and n r is the num- 
ber of correlated channels. The average center of mass 



o 

X 



Q 2 
1 





; (a) 






F = 0.05 




T = 0.0022 : 


L 


L = 60 



20 



40 



0.050 
0.045 : 
0.040 : 
Q 0.035 - 
0.030 
0.025 
0.020 






F = 0.05 
T = 0.0022, L = 60 



20 



40 



t 



FIG. 22: Center of mass diffusion constants (a) D yy and (b) 
D xx vs. time t for the smectic phase in the high drive limit, 
F = 0.05, T = 0.0022, for a system of size 60 x 60. That 
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FIG. 23: Center of mass diffusion constants (a) D yy and (b) 
D xx vs. time t for the liquid phase at, F — 0.10, T = 0.006, 
for a system of size 60 x 60. 



displacement in time t is then AX CI[L = (n r y/J L/ N c )Xt = 
{ n r IWJ L])\t, where we used N c = fL 2 is the total num- 
ber of charges. Because it is a Poisson process, the vari- 
ance of the number of avalanches is equal to the aver- 
age, and so the fluctuation about this average displace- 
ment is (AA" cm ) 2 = (n r /[yffL]) 2 \t. This yields the ra- 
tio D xx /v avcx = N c (AX cni ) 2 /2AX cm = n r V7V 2 - For 
/ = 1/25 and L = 60 we get D xx /v avcx — 6n r . At 
low T and F the correlation length £j_ can get large (see 
Fig. ll6|) : if several channels are correlated, the ratio 40 
can be attained. 

For the high drive case of Fig. 1221 we find the ratio 
D xx /v avcx ~ 0.00025. In this limit where rows of chan- 
nels move steadily forward one after the other, longitu- 
dinal fluctuations are greatly suppressed. Finally, for the 
liquid case of Fig.EU we find D xx /v avcx ~ 25. In the liq- 
uid, the system is structurally disordered and the motion 
of the charges is largely uncorrelated. Fluctuations about 
the center of mass motion are correspondingly enhanced. 



IV. RESULTS ON A SQUARE GRID 

We now consider the behavior of the driven Coulomb 
gas on a periodic square grid of sites. We consider only 
CTMC dynamics for the same charge density of / = 1/25 
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that was considered above for the triangular grid. We 
first consider behavior in the limit T — > 0. In Fig. 124b 
we show the equilibrium ground state configuration for 
F = 0. The charges occupy the sites of a 5 x 5 square sub 
lattice of the grid. The basis vectors of this sub lattice, 
ci = 3i— Ay and C2 = Ax+3y, are clearly not aligned with 
the grid basis vectors di = x and 02 = y, nor with the 
driving force F = Fx that we apply. This will produce 
some interesting effects. 

When F > F c ~ 0.06 the charges will start to move for- 
ward parallel to F, according to the order in which they 
most lower the system energy. In Fig. l24b we number 
the charges in the order in which they move in a par- 
ticular CTMC pass, and in Fig.l^b we give the change 
in interaction energy AH associated with each move, as 
was done for the triangular grid in Fig. [2 Charges move 
forward in the x direction in an order dictated by their 
position along the c = Ci + c 2 direction, as indicated by 
the arrow drawn in Fig. l24b . If one follows a path along 
the direction of c, using periodic boundary conditions, 
one finds that the path closes upon itself only after one 
has passed through all the charges in the ground state. 
Thus there is no row by row motion as there was for the 
case of the triangular grid, and hence no oscillation in 
AH as a function of simulation step. 



A. High Drive 

We now consider behavior at finite temperature and 
high drive. We simulate a 50 x 50 size system, starting 
from the F — ground state, at the values T = 0.004, 
F = 0.10. In Fig. 125b we show an intensity plot of the 
structure function S(k) at the initial stage of the simu- 
lation; our results are averaged over 1250 CTMC passes 
after an initial 2500 passes were discarded for equilibra- 
tion. We see peaks at wavevectors K corresponding to 
the ordered F = ground state. The peaks remain sharp 
in the k\ direction, but are somewhat smeared out in 
the transverse direction, suggesting a moving lattice 
with anisotropic translational correlations. If we simulate 
longer however, this moving ground state lattice under- 
goes a change of structure. In Fig.l25b we show 5(k) 
averaged over 5 x 10 6 passes, after discarding an initial 
5 x 10 6 passes. We see clearly a 6-fold orientational order 
in the position of the peaks, which are aligned with one of 
the diagonals of the square grid. Fig. 125b is reminiscent 
of the floating triangular lattice (algebraic translational 
order) that is seen in equilibrium simulation^ of more di- 
lute systems on a square grid, however without a finite 
size scaling analysis we cannot be certain of the nature of 
translational order in the system. Finally, however, if we 
simulate even longer, the structure changes yet again to 
a smectic phase with channels oriented parallel to F. In 
Fig. 123: we show 5(k) averaged over 2.5 x 10 7 passes, af- 
ter discarding an initial 3.75 x 10 7 passes. We see clearly 
the same smectic structure that we saw for the trian- 
gular grid in Fig.^Jl. The extremely long (compared to 
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FIG. 24: CTMC on a square grid with charge density / = 
1/25 at T — > and F > F c , with F parallel to the &i axis, 
(a) Ground state charge lattice for a 25 x 25 square grid. 
Numbers denote the locations of the charges in the ground 
state. The value of each number indicates the step on which 
that charge moves forward, (b) The change in interaction 
energy ATI at each step as charges move forward. (•) are for 
a 25 x 25 grid and correspond to the moves in (a); (o) and 
(a) are the beginnings of similar sequences for 50 x 50 and 
100 x 100 grids. 



the triangular grid) time it takes for the system to order 
into the smectic results because the initial ground state 
configuration is ordered with a set of reciprocal lattice 
vectors {K} that is not commensurate with that of the 
final state smectic. The system first requires a long time 
to disorder the initial state, and then another long time 
to reorder into the smectic. 

We have checked the above results by carrying out sim- 
ulations in which we start from an initial random config- 
uration of charges. After roughly 3.5 x 10 7 passes, we find 
that the system orders into the same smectic state as in 
Fie .l25b. In Fig.l25H. we show an intensity plot of the real 
space correlations in the smectic, C(mi, 777.2), obtained 
from the Fourier transform of the S(k) of Fig.l25b. We see 
that charges in the same channel (777,2 = 0) have a sharp 
periodic ordering. Correlations between channels show 
that the charges in neighboring channels are staggered; 
the peak in charge density in one channel aligns with the 
minimum in charge density of the neighboring channel, 
so as to form a local ordering that is more triangular than 
square. As one moves to channels further away, the corre- 
lations decrease and the peaks in C{m\, 7772) get smeared 
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FIG. 25: (a-c) Intensity plot of structure function S'(k) for a 
50 x 50 size system at T — 0.004, F = O.lSi , starting from the 
ground state of Fig. l24fa . (a) S(k) averaged over 3750 passes, 
after an initial 2500 passes of equilibration; (b) S(k) averaged 
over 5xl0 6 passes, after discarding an initial 5xl0 6 passes; (c) 
S(k) averaged over 2.5 x 10 7 passes, after discarding an initial 
3.75 x 10 7 passes, (d) Intensity plot of real space correlations 
C(mi,m2) corresponding to the smectic 5*(k) of (c). 



out. 

We now check that the smectic phase of Fig.l25b has 
the same scaling behavior with system size that was 
found for the triangular grid. For larger systems with 
L = 100 — 200, it is not possible to simulate for the very 
long times (~ 10 7 passes) that are needed to order into 
the smectic from either the ground state or a random 
initial state. We therefore start with an initial configu- 
ration that is a periodic repetition of the smectic state 
for L — 50, and simulate for only relatively short times. 
For L = 100, 150, and 200, we use 10000, 2000, and 
2000 passes. Our results are shown in Fig.[2BI where we 
plot the profiles of S(k) along different paths through 
the first Brillouin zone. Fig.l2rJb shows that the speaks 
are as sharply confined to the values k\ = 1/5,2/5 as 
was found for the triangular grid. Fig. 120b shows that 
the peaks at k% = scale ~ L 2 , indicating long range 
smectic order. Figs.l2l)b.d. show that S(kx,k2) for fixed 
k% = 1/5,2/5 scales rougly ~ L. Note that the scaling 
collapse in Fig.l26b is not quite as nice as the correspond- 
ing Fig.0; for the triangular grid. Plotting the peak value 
S(Kii) versus L, similar to what was done in Fig. El gives 
S(Kii) ~ L s with s w 1.17. We believe that this value 
s > 1, rather than being a signature of stronger correla- 
tions between smectic channels, may just reflect the per- 
sistence of correlations introduced by our initial periodic 
configuration, which have not yet completely washed out 
over our relatively short runs. 

In Fig.|^|we plot the transverse and longitudinal cor- 
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FIG. 26: Profiles of S(k) in various directions, for different 
system sizes L, for the smectic phase at F = 0.10, T = 0.004 
on a square grid, (a) 5*(k) vs. ki for fixed ki = 1/5; (b) 
5(k)//L 2 vs. k 2 for fixed ki = 0; (c) S{k)a /L vs. k 2 for 
fixed fei = 1/5; (d) S(k)a /L vs. k 2 for fixed ki = 2/5. Note 
the logarithmic scale in (a) and (b). 
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FIG. 27: Transverse and longitudinal correlation functions at 
F — 0.1, T = 0.004, for system sizes L x L on a square grid, 
(a) \C(ki — 1/5, at values 7712 = 5n, n integer, for differ- 
ent sizes L; solid lines are fits to the periodic exponential of 
Eq. 13511 . (b) C(mi, m2 = 0) at values mi = 5n, n integer, for 
different sizes L; solid lines are fits to the periodic exponential 
of Eq. 



relation functions, obtained from the appropriate Fourier 
transform of S(k). C(fci = 1/5, mi) in Fig.l27b shows ex- 
ponentially decaying transverse correlations, with a cor- 
relation length £j_ ^5 — 7, comparable to that found 
for the triangular grid at the same parameter values. 
We believe that the slight increase from £j_ ~ 5 to 7 
as L increases from 50 to 200 reflects the correlations 
introduced by our initial configuration, as already com- 
mented on in connection with Fig. 126b. Note that the 
peak values of C{k\ — 1/5, m^) oscillate in sign for suc- 
cessive values of 7712 = 5n, n integer, due to the ir phase 
shift from channel to channel that is apparent in the real 
space correlations shown in Fig. l25t h hence we have plot- 
ted \C(ki = 1/5, m 2 )| in FigJSBk. In Fig.E7b we plot the 
longitudinal correlation C{m\,m2 = 0). The solid lines 
are fits to a periodic exponential, and give a common 
value of £|| ~ 174 for all sizes L. Thus we have a finite 
correlation length, but £|| ~ L. We conclude that the 
driven steady state for low T and large F on the square 
grid is a smectic that is qualitatively the same as what 
was found for the triangular grid. 
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FIG. 28: 6-fold orientational order parameter | «3>6 ] vs. simu- 
lation clock time t for T = 0.004, F = 0.04 and L x L system 
of size L — 75. Regions denoted "L" are in a moving liquid 
state; regions denoted "SI" and "S2" are in a moving solid 
state. Coexistence of the two states indicates that the system 
is at the melting transition. 



B. Low Drive 

We now consider behavior at low drive, simulating at 
T = 0.004, F = 0.04 for anLxL system of size L = 75. 
We will find that these parameters place the system right 
at the melting transition. We start from an initial ran- 
dom configuration and run 2.5 x 10 4 passes to equili- 
brate, followed by 2.5 x 10 7 passes to compute averages. 
In Fig-EHl we plot the instantaneous absolute value of 
the 6-fold orientational order parameter |$g| versus the 
simulation clock time t. We see that the system makes 
sharp jumps between states of lower and higher values of 
I $6 1 and conclude that these are the coexisting liquid and 
ordered phases at the first order melting transition. In 
Fig.l29b we show an intensity plot of the structure func- 
tion S(k) averaged over only the liquid states labeled "L" 
in Fig. ESI We see a liquid like S(k), but with a striking 
6-fold modulation of intensity in the diffuse peaks, corre- 
sponding to the relatively large values of |$g| ~ 0.4 seen 
in Fig-EHl In Fig-I29b we show 5(k) averaged over the 
ordered states labeled "S2" in Fig.[2Hl We see periodic 
sharp peaks suggesting a moving solid state. S(k) for 
the states labeled "SI" in Fig.[2Hlis identical to that of 
Fig.l2l3b. except reflected about the k 2 axis. In Figs.l23b.d 
we show intensity plots of the corresponding real space 
correlations C(r). 

We consider first the liquid state. In Fig.l30b we plot 
S(k±,k2 — 0) versus ki for several different L x L sys- 
tem sizes. Except for the maximum of the first peak, 
we see essentially no finite size effect. The height of the 
first peak is different for the different L, however there 
is no systematic variation with L; we believe that these 
differences are just statistical fluctuations. We conclude 
that this state is a liquid with only short range transla- 
tional order. Figs .12*51 andl2"9"k however suggest that the 
liquid may possess finite orientational order. In Fig. l30b 
we therefore plot |($6)| versus L for the liquid state. We 
see that |($6)| is roughly independent of L, confirming 
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FIG. 31: Schematic showing how particle clusters with local 
6-fold orientational order (shaded circles connected by thin 
lines) may align with the rotational symmetry of an external 
potential or grid (thick lines), (a) lock in of the cluster to a 
6-fold rotationally invariant triangular grid; (b) lock in of the 
cluster to the horizontal axis of a 4-fold rotationally invariant 
square grid; and (c) lock in to the vertical axis of a square 
grid. 



FIG. 29: Intensity plots of S{k) at T = 0.004, F = 0.04, for an 
LxL system of size L — 75 averaged over (a) the liquid states 
labeled "L" and (b) the solid states labeled "S2" of Fig-EHl 
Instensity plots of the corresponding real space correlations 
C(r) for (c) the liquid states "L" and (d) the solid states "S2". 
The applied force F is in the horizontal direction. 
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FIG. 30: The liquid state at T = 0.004, F = 0.04 for different 
LxL system sizes, (a) S(ki,k2 = 0) vs. ki; (b) 6-fold 
orientational order parameter |($6)| vs. L; bars denote the 
standard deviation of the distribution of |$6|- 



that the liquid has long range hexatic orientational order. 

The long ranged hexatic liquid that we find in Fig. 129b 
is reminiscent of the long ranged hexatic liquid found 
for the triangular grid at low temperatures, as shown in 
Fig. 1181 There are, however, some important differences. 
A liquid in a continuum will always have local 6-fold ori- 
entational order. However, due to the short ranged trans- 
lational correlations, the phase of the local complex ori- 
entational order parameter will vary with position. For 
a normal liquid, this causes correlations of the 6-fold ori- 
entational order parameter to decay exponentially with 
distance. According to the theory of melting in two di- 
mensions by Halperin and Nelson, and by Youngi^, there 
may also be an algebraically ordered hexatic liquid be- 
tween the solid and normal liquid phases, in which corre- 
lations of the 6-fold orientational order parameter decay 



algebraically. When the system sits on an external pe- 
riodic potential, however, the local 6-fold orientational 
order parameter can lock onto the symmetry directions 
of the external potential, which therefore serves as an 
ordering field for orientational order. For a triangular 
grid, the local 6-fold order of the particles locks onto the 
6-fold rotational symmetry of the grid, as illustrated in 
Fig.lblb. The result is long range 6-fold orientational or- 
der, with a finite ($6)- For a square grid, the local 6-fold 
order of the particles may lock onto either the vertical 
or the horizontal directions of the grid, as illustrated in 
Figs. liSlb .c. For a liquid in equilibrium, both of these ori- 
entations will occur in equal numbers on average. Since 
the two orientations are related by a it/ 2 rotation, the 
relative phase of the 6-fold orientational order parameter 
for the two cases is exp(z67r/2) — — 1, and adding them in 
equal numbers causes (^e) to vanish. A square grid in- 
duces no 6-fold orientational order in equilibrium. For a 
liquid in a driven non-equilibrium steady state, however, 
the direction of the driving force F breaks the symmetry 
between the vertical and horizontal directions, and can 
cause one to be favored over the other. A driving force 
therefore can lead to a finite ($g) and long range 6-fold 
orientational order on the square grid. From the plot of 
the real space correlation C(r) shown in Fig.l29b. we see 
that the system locks onto the vertical direction, as in 
Fig. l31b . The resulting structure function S(k), shown 
in Fig.l29b has a set of 6 peaks about the origin, which 
are oriented so that one pair of the peaks align with the 
direction parallel to the applied driving force F. This 
is in contrast to case for the triangular grid, shown in 
Fig-El where the peaks are oriented so that one pair 
of the peaks align with the direction transverse to the 
applied force. 

We now consider the ordered moving state, labeled 
"S2" in Fig.EHl The structure function 5(k), and the 
real space correlations C(r) are shown in Figs. 129b and 
1231 respectively. Note that the periodic peaks in 5(k) 
for this ordered state do not have the same symmetry 
as that of the equilibrium ground state. The latter (see 
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Fig.l25b) consists of a square array of Bragg peaks, while 
in Fig.l29b the peaks are distorted into a more triangular 
structure. From either the location of the peaks in 5(k), 
or more easily from a direct inspection of the real space 
correlation C(r), we see that this state consists of peri- 
odic channels of charges oriented parallel to the applied 
force F in the a\ = x direction. Within each channel 
the charges are ordered with an average separation of 
8^ grid spacing, while the channels themselves are sep- 
arated from each other by 3 grid spacings. The nearest 
neighbors to a given charge are located in its two nearest 
neighboring channels, rather than within the same chan- 
nel, reflecting a similar orientation of hexatic order as in 
the liquid. This can be compared with the smectic state 
at high drive, shown in Fig.l25H. This smectic has chan- 
nels in which charges are separated by 5 spaces, while 
the channels themselves are separated by 5 spaces; the 
nearest neighbors to a given charge are located within 
the same channel. In contrast, the equilibrium ground 
state of Fig.lMt can be thought of as channels in which 
charges are separated by 25 spaces, while the channels 
themselves are separated by 1 space; nearest neighbor 
charges lie in the next-next-nearest neigbhoring channels. 
The structure in terms of channels, as described above, is 
determined by the strength of the correlations between 
the channels. When channels arc more strongly corre- 
lated, it can be energetically favorable to have the chan- 
nels spaced more closely together, with a correspondingly 
larger distance between charges within a given channel; 
the stronger correlations between the channels will keep 
charges within neighboring channels from approaching 
each other too closely. The equilibrum ground state rep- 
resents the extreme limit of long range correlations be- 
tween channels. The high drive smectic represents the 
opposite limit where correlations between channels are 
very short ranged and it becomes favorable to keep the 
spacing between channels at the same distance as the 
spacing of charges within a channel. The moving state 
of Fig.l29tl can thus be thought of as having a struc- 
ture, and presumably channel correlations, intermediate 
between these two limits. 

The strong correlations between channels, as discussed 
above, are clearly evident in the plots of 5(k) and C(r) 
in Figs.l2*% and [29H . All peaks in S(k) appear sharp 
in both the longitudinal and transverse directions; real 
space correlations C(r) appear to extend the entire length 
of the system in both longitudinal and transverse direc- 
tions. This suggests a moving solid rather than a smec- 
tic. To investigate this further we consider in more detail 
the peaks in the structure function 5(k). In Fig.l32b we 
plot profiles of S(ki, fe) versus k2, showing the peaks at 
fci = and fci = 3/25. For k\ = the peaks appear 
as sharp 5-function peaks upon a smooth background, 
similar to what was seen in Fig.l26b for the smectic at 
large drive, indicating the periodicity of the channels in 
the direction transverse to the driving force F. At finite 
fci = 3/25 the peaks are much sharper than the corre- 
sponding finite fci peaks for the smectic at large drive 
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FIG. 32: Moving ordered phase of Fig.l^b. for T = 0.004, 
F = 0.04, and L = 75. (a) Profiles of S(ki,k2) versus k 2 
showing the peaks at fci = and fei = 3/25; note the loga- 
rithmic scale, (b) Heights of the dominant peaks in S(ki, k 2 ) 
versus k 2 ; the different curves represent different values of fci 
and the solid lines are fits to a Gaussian. 
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FIG. 33: (a) Moving ordered phase of Fig.l2"5b. for T = 0.004, 
F = 0.04, and L = 75. Values Sm(ki,k 2 = 0) of fits from 
Fig.CHt vs. fci. Solid line is a fit of the data to a Gaussian, 
(b) Moving ordered phase for T = 0.003, F = 0.04, and sizes 
L = 75, 150, and 225. Values of S Rt (ki,k 2 = 0) vs. Lk\. 
Solid line is a fit of the small k\ data to a straight line. 



in Fig.l26b: in the present case the peaks drop by three 
orders of magnitude from maximum to minimum (note 
the logarithmic scale) and have a half width of about 
Afc ~ 0.007. Such sharp peaks suggest the possible pres- 
ence of long ranged or algebraic correlations between the 
channels. 

In Fig.l32b we plot only the heights of the dominant 
peaks in S{k\,k<i) versus k%, for the different values 
of fci. We see that at fixed fci, there is only a very 
small variation of the peak heights with fc2 ; however the 
dependence on fci is considerable. Fitting the points 
for each value of fci to a simple Gaussian (the solid 
lines in Fig.l52b1 we plot the resulting Sat (fci, &2 = 0) 
versus fci in Fig. l33b . where another simple Gaussian, 
Sfit(fci, 0) = A c exp(— afcf), gives an excellent fit (the 
solid line in Fig. 133b) . Such a Gaussian shape for the 
peak heights is consistent with a Debye- Waller-like be- 
havior for thermal fluctuations of a solid. 

However to investigate more precisely the nature of 
translational correlations, we need to investigate the de- 
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pendence of the peak heights on system size L. We have 
not, however, been able to do this at T = 0.004; when, 
for larger system sizes L, we start the system off in an 
initial disordered state, we were unable to see a similar 
transition to an ordered state as was found in Fig. 1281 for 
L = 75. We assume that this is either because the melt- 
ing temperature becomes somewhat lower for larger L (as 
was seen for the triangular grid) , or perhaps because the 
free energy barrier between the liquid and ordered states 
increases with L and we have not run sufficiently long to 
have a thermal excitation over this barrier. 

We choose, therefore, to investigate the finite size be- 
havior at the lower temperature T — 0.003, taking as an 
initial state an appropriate cut out of, or periodic exten- 
sion of, the ordered state we found for T — 0.004. For 
sizes L — 50 and L = 100, when we started the sys- 
tem in such an initial state, we found that the system 
quickly melted to a liquid. We believe that this is be- 
cause these values of L are not commensurate with the 
spacing of 3 grid spaces between channels required by 
this ordered structure. Ordered systems of size L = 75, 
150 and 225, however, remained stable. Proceeding simi- 
larly to Fig.l52b. we examin the peak heights of S(ki, k%) 
versus ki for the various k\. At this lower temperature 
T = 0.003, the variation with fc 2 is even smaller than 
that seen in Fig.l3*2*b at T = 0.004. Fitting to a Gaussian, 
we determine the values of Sfn(ki, = 0) and fit these 
to a Gaussian, Sst(kx,Q) = N c exp(— a(L)k\), where 
N c = fL 2 . We find the surprising result that a(L) ~ L. 
To show this, we plot in Fig.|2St> Sst(ki,0)/N c , on a log 
scale, versus Lk\. We see that the data at small k\ give 
an excellent collapse to a straight line. This exponential 
decrease in S(K)/N C with increasing L suggests that the 
observed moving solid may not persist as a stable state 
in the large L limit. 

We can also examine the translational order from the 
perspective of the real space correlations. In Fig. 134b we 
plot the longitudinal correlations, C(toi,TO2 = 0) versus 
mi, for system sizes L = 75, 150 and 225. We clearly see 
that there are three charges for every 25 grid spacings. 
No finite size effect is seen. In Fig. 134b we plot the abso- 
lute value of the complex correlation C(k 1 = 3/25, 7712) 
versus mi , showing values for only every third grid spac- 
ing, 7TI2 = 3n, n integer. Here we find a pronounced finite 
size effect, with the correlation decaying to lower values 
as L increases. However a periodic exponential (as used 
for example in Fig.l27b1 does not give a particularly good 
fit, and we do not have enough sizes L to try any sys- 
tematic scaling fit. While the results of Figs.l3*5b and!34b 
thus suggest that long range solid order may not persist 
as L increases, larger sizes will be needed to clarify the 
true large L behavior. 

The structure that we have found in Fig.l29b.d for the 
ordered moving state at low drive has neither the com- 
mensurability with respect to the underlying grid of the 
equilibrium ground state, nor the large drive smectic. 
One can speculate that at other values of F and T, in 
this low temperature ordered region, yet other commen- 
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FIG. 34: Real space correlationf for the moving ordered phase 
at T = 0.003, F = 0.04, and L = 75, 150, 225. (a) Longi- 
tudinal correlation C(mi,m2 = 0) vs. mi. (b) Transverse 
correlation \C(ki = 3/25, ma) | vs. m.2, for 7712 = 3n, n inte- 
ger. 



surabilities may be found. Exploring the complete phase 
diagram of the driven lattice Coulomb gas on the square 
grid may therefore prove to be considerably more chal- 
lenging than was for the case of the triangular grid, and 
we leave this for future investigations. 



C. Dynamics 

We now consider some of the dynamic properties for 
the driven Coulomb gas on the square grid. We will not 
attempt a detailed calculation of diffusion constants, as 
we did for the triangular grid, however we will still be 
able to make some interesting observations by looking 
at average velocity and center of mass motion. We first 
consider the case of the high drive smectic, F = 0.10, 
T = 0.004, considered in section llV Al In TablelTTlwe give 
the values for the average center of mass velocity parallel 
to the driving force, v avox , for various system sizes Lx L. 
Similar to our results for the triangular grid (see Table HJ) 
we find i> avc x ~ L scales proportional to the length of the 
system in the direction of the applied force, in agreement 
with the discussion at the end of section lll B 31 Inspection 
of the center of mass motion as a function of time clearly 
shows no transverse diffusion, indicating that the smectic 
is transversely pinned, just as we found for the triangular 
grid. 



TABLE II: Average velocity u avcl in the high drive smectic 
state for various system sizes L x L at F = 0.10, T = 0.004, 
on the square grid. 



L 


100 


150 


200 


^avc x 


1082 


1617 


2065 



Next we consider the case of low drive, F = 0.04, con- 
sidered in section lTV Bl We consider first the case at melt- 
ing, T — 0.004 and L — 75, where the system is making 
transitions between the liquid and a more ordered state, 
as shown in Fig. 1281 In Fig. 135b we plot the component 
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FIG. 35: Center of mass displacement for F — 0.04, T — 
0.004, and system size L = 75, at the melting transition: (a) 
motion parallel to F, X cm vs. time t; (b) motion transverse to 
F, Ycm vs t. Light lines correspond to times when the system 
is in the liquid state. Heavy lines correspond to times when 
the system is the ordered state (see Fig. 1281 . 



of the instantaneous center of mass displacement paral- 
lel to the driving force, X cm , versus the simulation clock 
time t. Light lines denote times in which the system is in 
the liquid state, while heavy lines denote times when the 
system in the ordered state (compare with Fig.l2*81). That 
the lines in each region of time are perfectly straight in- 
dicates a constant average velocity u avc x m each region. 
Note, however, that the velocity in the ordered state is 
slightly larger than that in the liquid state: in the liquid, 
v avC x = 10.05, while in the ordered state, v avcx = 11.50, 
or 14% larger. In Fig. 135b we plot the transverse compo- 
nent of the center of mass displacement, Y cm , versus the 
simulation clock time t. In the liquid, Y c shows the noisy 
fluctuations characteristic of diffusion; the observed bias 
towards increasing values of Y cm we believe is just a sta- 
tistical fluctuation. In the ordered phase, however, Y cm 
stays essentially constant indicating that the system is 
transversely pinned. 

We now consider the ordered phase at the lower tem- 
perature T = 0.003. In Table ITTT1 we give the values of 
!J avcl for systems of different sizes L x L. We consider 
only the values L = 75, 150 and 225 that result in an or- 
dered moving state. In contrast to the high drive smectic 
(see Table ITTf) we now find that v avex is independent of L. 
This at first seems paradoxical, since the ordered state at 
low drive is more strongly correlated than the smectic at 
high drive. However it may be that the incommeasura- 
bility of the ordered state in the parallel direction (where 
the average spacing between charges is 8^ grid spacings) 
is sufficient to remove the energy barriers responsible for 
the avalanche effects (see section III B 3D that give rise to 
the «avc z ~ L\ dependence in the high drive smectic. 
Finally, an examination of the transverse displacement, 
similar to that of Fig.l35b. shows that the ordered state 
is transversely pinned. 



TABLE III: Average velocity « avcl in the ordered phase at 
various system sizes L x L for F = 0.04, T = 0.003, on the 
square grid. 



L 


75 


150 


225 


^avc x 


36.3 


36.3 


36.1 



V. DISCUSSION AND CONCLUSIONS 

In this work we have applied lattice gas dynamics to 
model the non-equilibrium steady states of a driven dif- 
fusive system, the 2D classical lattice Coulomb gas in 
a uniform applied force. We have considered two dif- 
ferent dynamic algorithms, and have found that they re- 
sult in qualitatively different phase diagrams, contrary to 
naive expectations. We have shown that the commonly 
used driven diffusive Metropolis Monte Carlo algorithm 
(DDMMC) results in a structurally disordered moving 
steady state over most of the phase diagram. We have 
argued that this is due to unphysical intrinsic random- 
ness in the algorithm that remains even as T — > 0. 

We have then applied continuous time Monte Carlo 
(CTMC) to the driven diffusive problem and found it 
to produce smectic and, for the square grid, possibly 
more strongly correlated steady states at low temper- 
atures. We have shown (Appendix A) that CTMC is 
a natural discretization of continuum Langevin dynam- 
ics. We have argued (section III B 2|) that in general it 
gives a physically correct dynamics when the grid sites 
are regarded as the minima of an external one body po- 
tential U(r), and the energy barriers Uq of this potential 
remain larger than the energy change AE in hopping 
between neighboring minima, so that motion is by ther- 
mal activation of one particle at a time over the poten- 
tial barriers. It remains unclear whether or not CTMC 
will be qualitatively correct in the very large drive limit, 
AE > Uq, when the applied force overcomes the pinning 
force of the potential, and the minima of the correspond- 
ing washboard potential, U(r) — F • r, vanish. In such a 
case, for a spatially uniform system in an initially ordered 
state, each charge will experience an equal net force for- 
ward from the washboard potential, and one would ex- 
pect at low temperatures that the charges would move 
coherently together. The CTMC algorithm, which only 
moves a single charge at a time, breaks this spatial uni- 
formity and might introduce unphysical effects. For the 
case of a system with quenched randomness, however, 
the random pinning already breaks spatially uniformity, 
and the forces on the charges will in general be differ- 
ent. In such a case, the single particle moves of the 
CTMC algorithm may not be as unphysical. This very 
large drive limit for the case of random pinning has been 
the subject of numerous recent theoreticali£*i2i22i2i and 
numerical 17 i 21 i 22 i 23 i 24 i 25 i 26 i 27 works. 

For CTMC we have shown (section [HI Bp that diverg- 
ing correlation lengths as T — > can give rise to subtle 
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finite size effects that can be difficult to detect with the 
usual finite size scaling methods applied to the peaks of 
the structure function, S(K), and we have argued that 
the smectic state that we find for finite size systems will 
become unstable to a liquid on sufficiently large length 
scales. However, since the relevant correlation lengths 
diverge as T — ► 0, the smectic will be the stable steady 
state in any finite system, at sufficiently low temperature. 
We have also shown that, on a square grid, long range 
hexatic orientational order develops in the moving steady 
state liquid, and that this is a purely non-equilibrium ef- 
fect. 

The one component 2D lattice Coulomb gas serves as 
a model for logarithmically interacting point vortices in a 
2D superconducting network, or a superconducting film 
with a periodic potential. Driven vortices in a 2D peri- 
odic potential at finite temperature have been simulated 
by several others using continuum dynamics. The molec- 
ular dynamic simulations of Reichhardt and Zimany: 28 
and of CarneiroSi used square periodic pins embedded 
in a flat continuum, with a number of vortices equal to, 
or greater than, the number of pins. We do not expect 
that such models, in which a sizable fraction of the vor- 
tices spend most of their time in the flat space between 
the pins, will be well described by our dilute density of 
charges on a discrete grid, where all charges spend most 
of their time at the potential minima. 

Much closer to our model is that of Marconi and 
Dommguea22i2i who simulate the dynamics of a square 
array of Josephson junctions using resistively-shunted- 
junction (RSJ) dynamics applied to a 2D XY model. 
They study a vortex density per unit cell of the array of 
/ = 1/25, the same density as used in our present work. 
The phase diagram which they report has some quali- 
tative similarity to our own phase diagram of Fig.Qb, 
with an ordered, transversely pinned, moving state at 
low temperatures^ 2 -. However, in contrast to either the 
smectic we find at high drives (see Fig.l25b) or the more 
ordered state we find at low drives (see Fig.l29b). they 
find a moving vortex lattice where S(k) has peaks at 
the same reciprocal lattice vectors K as the equilibrium 
ground state. From a finite size scaling analysis of 5(K) 
using L — 50, 100, and 150, they conclude that their state 
is a vortex lattice with anisotropic algebraically decaying 
translational correlations. 

To understand possible reasons for the difference be- 
tween their results and ours, we first consider the relevant 
parameters of their model. For their cosine Josephson 
junction model, the effective one body potential 33 that 
the array structure introduces for vortex motion has an 
energy barrier Uo — 0.12, in units where the Joseph- 
son coupling energy is Jo = 1. Many of Marconi and 
Dominguez's results are in the limit where the force F 
(i.e. the applied current in the Josephson array model) 
satisfies F > Uo. This is the case where the minima 
in the washboard potential parallel to the driving force 
have vanished, and where we have argued that our lat- 
tice gas dynamics might not apply. However, even for 



the case F < Uq, the two models may be in different 
parameter regimes. Our lattice gas dynamics implic- 
itly assumes that the energy barrier Uq is larger than 
all other energy scales. For the Josephson array of Mar- 
coni and Dominguez, however, a direct calculation using 
the XY model shows that the energy to move a single 
vortex forward one grid space from its ground state posi- 
tion is AE\ ~ 0.34, substantially bigger than the barrier 
Uq — 0.12. Our simulations are therefore in the limit of 
a much stronger pinning potential. 

In spite of these parameter differences, we can still 
make some observations. First we note that Marconi 
and Dominguez always begin their simulations from the 
equilibrium ground state (or states evolved from it); they 
are unable to cool the system from a liquid and find the 
ordered state, hence there is no independent check that 
the state they find is the true stable steady state. Next, 
we note that because their simulations use a continuum 
dynamics, they are unable to simulate for the very long 
times that are possible using our lattice gas dynamics. 
As a measure of the effective simulation time, we can 
compute the total displacement Ai? cm of the vortex cen- 
ter of mass parallel to the applied force over the total 
time of the simulation. For the Josephson array, if V 
is the average measured voltage drop per junction par- 
allel to the applied current, Jo the critical current of a 
single junction, Rn the normal shunt resistence, / the 
vortex density, and tj = h/ {2eR^Io) the time constant, 
then AR cm = (V/IoR N )(At/Tj)(N t /2nf), where At is 
the time integration step of the simulation, and N t is the 
number of such steps. Using Marconi and Dominguez's 
values 3 ^ of At/rj = 0.1, / = 1/25, N t = 10 5 , and typical 
values^ of V/IqRn from their Fig. 5, we find for their 
simulations that Ai? cm ~ 1.2 x 10 3 grid spacings or less. 
In contrast, our simulations which lead to Fig.l25h have a 
total simulation time corresponding to Ai? cm ~ 6 x 10' 
grid spacings, more than 4 orders of magnitude larger. To 
make a better comparison with Marconi and Dominguez, 
we note that our results of Fig. l25h ,. starting from the 
equilibrium ground state, correspond to a total center of 
mass displacement of only Ai? cm ~ 4 x 10 3 grid spacings, 
similar to that of Marconi and Dominguez. The state we 
find in Fig. l25h has peaks in 5(k) at the same K as the 
equilibrium ground state, moreover the anisotropies of 
this state are the same as for the state found by Marconi 
and Dominguez; the peaks develop a finite width in the 
direction transverse to the direction of motion (this fea- 
ture is visible in Fig.|53^), and the heights of the peaks 
decrease as k varies in the direction of motion. In our 
case the variation in peak heights is only a 20% reduc- 
tion from largest to smallest, whereas for Marconi and 
Dominguez it is a larger 75%, nevertheless the behavior 
is qualitatively similar. It thus may be that the sim- 
ulations of Marconi and Dominguez have not run long 
enough to observe the true long time steady state of the 
system. 

Finally, we comment on one additional issue that is 
related to our ability, using lattice gas dynamics, to sim- 
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ulate to much longer times that can be achieved with 
continuum methods. It is interesting to note in Fig. 1231 
that the longitudinal diffusion constant D xx in the liq- 
uid approaches its long time limit on a much longer time 
scale than does the transverse diffusion constant D yy . In 
recent continuum Langevin simulations 34 of driven vor- 
tices in a disordered 2D superconductor, similar diffusion 
in the vortex liquid phase was computed. Although it was 
observed that the transverse diffusion constant saturated 
to a finite value at long times, the longitudinal diffusion 
constant was found not to saturate, but rather to grow 
proportional to t. Rather than reflecting super-diffusive 
behavior in the longitudinal direction^*, this result might 
simply be a failure to simulate to long enough times to 
see the longitudinal diffusion constant saturate. 
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Next we symmetrize the Fokker-Planck equation by 
making the transformation, 



^({r,}) = e^>]/^P( {l .}) 



(44) 



Substituting the above into the Fokker-Planck equation 
(|43|l gives the imaginary time Schrodinger equation, 



(hp 

at 



d 2 ip 
dr? 



where, 



V in = 



1 dU 
2T^~ 



Vi a 1p 



1 d 2 U 
2Tdr!~ 



(45) 



(46) 



If we now discretize the coordinates, so that the r; are 
confined to the sites of a periodic lattice, the natural way 
to discretize Eq. (|45|l is to replace the second derivative 
of ip with its lattice equivalent, 



dt 4^ 



'A? - Vi] i> 



(47) 



Appendix A 

In this appendix we demonstrate that the transition 
rates of Eq. i|13fl correctly describe diffusive Langevin mo- 
tion in the limit that the energy change in one move sat- 
isfies AU <C T. Our derivation follows one given earlie^ 
for a single degree of freedom, and extends it to the case 
of many degrees of freedom. 

The Langevin equation of motion for diffusively mov- 
ing particles in a uniform driving force F can be written 
as, 



dr. 



dU 



dt - D ^- + ^> (39) 
where Ti a is the a component of the position of particle 



[/[{r,}] ^ H[{r,}] - F • ^> 



(40) 



with TL the Hamiltonian giving the internal interactions 
between the particles, and rji a is the a component of the 
thermally fluctuating force acting on particle i. In order 
that the system reaches equilibrium in the absence of the 
force F, the thermal force is taken to have correlations, 







2DT5ij6 a p5(t - t') 



(41) 
(42) 



The corresponding Fokker-Planck equation that de- 
scribes the probability P({ri}) for the system to be at 
coordinates {r^} is then given by, 



dt ^ 



i) 



P 



dU 



dr ia V dr. 



T 



d 2 P 



drl 



(43) 



where Af is the discrete Laplacian for the lattice with re- 
spect to coordinate r i; and Vi the appropriate discretiza- 
tion of J2 a Via > as w i n ^ e explained below. For a lattice 
with nearest neighbors given by the vectors {ia^}, the 
discrete Laplacian acting on a scalar function /(r) is de- 
fined by, 

A 2 /(r) = c [/(r + S M ) - 2/(r) + /(r - a,)] , (48) 



with c an appropriate geometrical constant to give the 
correct continuum limit. 

If we denote the state of the system with particles at 
positions {ri} as s, then we can write the above Eq. l|4T|) 
in a matrix form, 



dips 
dt 



(49) 



where the matrix M has elements, 

M ss = -DT[zc+Vi\ (50) 
M ss , = cDT, whens' = {ri±a M ,rj} (51) 
M ss i = 0, otherwise , (52) 

where z is the number of nearest neighbor sites. By the 
notation in Eqs. H51J) and l|52(l we mean that the only 
non-zero off-diagonal elements of M are those connect- 
ing states s and s' in which only a single particle at has 
moved to a nearest neighbor position ri±a M , and all other 
particles have remained unchanged. This is our first re- 
sult: the natural discretization of continuum Langevin 
dynamics to a lattice gas dynamics involves single parti- 
cle moves only. 
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To see what are the correct single particle hopping 
rates for our lattice gas dynamics, as well as to see what 
is the correct discretized form for the Vi of Eq. (|47[) . con- 
sider now the Master Equation for our lattice gas dy- 
namics. If s = {ri} is the state of the system, then the 
probability P s to be in state s is determined by, 

dP 



where W s > s is the rate to make a transition from state s 
to state s'. We therefore have, 



M, 



M a8 , = W ss >, s^s' 



(54) 
(55) 



We next apply the same transformation as in Eq. (|44|l , 

to get, 



dips 

at 



Y,e u ° /2T M ss> e 



U S ,/2T 



1ps 



Comparing with Eq. (|49J) we get, 



(56) 



(57) 



and from Eqs. (|51J) and (|55(l we then get for the off- 
diagonal elements of M, 



M„, = e^- u ^y 2T W^ = cDT 



(58) 



when the state s' differs from the state s by only a sin- 
gle particle that has moved to a nearest neighbor site, i.e. 
r\; — > Ti±a^, with all other Yj kept constant; all other off- 
diagonal terms vanish. The above result then determines 
the transition rates that are needed for the discrete Mas- 
ter Equation to model the continuum Langevin equation, 



W ss , = cDTe-^~ u ^/ 2T 



(59) 



Thus we arrive at the rates of Eq. H13[) that define our 
CTMC algorithm. 



Note that the rates of Eq. (|59|l satisfy a local detailed 
balance, 



Wss' _ (U 3 -U 3 ,)/T 
Ws's 



(60) 



Having determined the rates W ss ' , we can now deter- 
mine the diagonal part of M and hence the Vi of Eq. 147|l . 
Defining AUi±^ as the change in U when a single particle 
moves ri — ► Ti ± dfj_, i.e., 



AU i±fl = U[{vi ± a^rj})] - U^Vj}] , 



(61) 



then Eqs. (|54ll and l|57|) give for the diagonal elements of 



Ms 



M ss = M ss = - E Ws's 
s' 

= -cDT^~ {U ''~ U3)/2T 

s' 

= - cDT H [ 



(62) 
(63) 
.(64) 



If one now expands Eq. (gU) to order (AU/2T) 2 , and then 
uses Eq. gEJ that cf^^AUi+f, + AL/j_ M ) = A 2 U, and 
compares to Eq. (|50")l . one concludes that, 



Vi 



Aft/ 
2T 



E 



c / AUi+,\ 2 , c {AU t - N ' 



2 V 2T 



2 V 2T 



This is just the natural symmetric discretization of Vi = 
E a v ™ with V a given by Eq. ®. 

We have thus shown that the CTMC dynamics, with 
rates as in Eq. (|13fl . is the natural discretization of over- 
damped Langevin dynamic in the continuum, and that 
CTMC becomes a very good approximation for the con- 
tinuum dynamics in the limit that the energy changes for 
single particle moves, AUi a of Eq. H61fl . satisfy AUi a -€i 
2T. 
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